# Nonlinear Grey Prediction Model with Convolution Integral NGMC(

Authors: Zheng-Xin Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The grey prediction model with convolution integral GMC (1,n) is a multiple grey model with exact solutions. To further improve prediction accuracy and describe better the relationship between cause and effect, we introduce nonlinear parameters into GMC (1, n) model and additionally apply a convolution integral to produce an improved forecasting model here designated as NGMC (1, n). The model solving process applied the least-squares method to evaluate the structure parameters of the model: convolution was used to obtain an exact solution with this improved grey model. The nonlinear optimisation took the parameters as the decision variables with the objective of minimising forecasting errors. The GMC (1, 2) and NGMC (1, 2) models were used to predict China’s industrial SO2 emissions from the basis of the economic output level as the influencing factor. Results indicated that NGMC (1, 2) can effectively describe the nonlinear relationship between China’s economic output and SO2 emissions with an improved accuracy over current GMC (1, 2) models.


---

## Body

## 1. Introduction

Due to the complexity of both the internal and external environmental factors in any such system, the behavioural data are often sparse. Statistical analyses can effectively address the problem where datasets are large; for small sample sizes (n
<
15), solution becomes difficult. Grey system theory, pioneered by Deng in 1982 [1], is an uncertainty theory dealing specifically with the analysis, modelling, prediction, and control of information-poor systems. Grey theory considers that although the objective system appears complex, with sparse data, it always has an innate overall governing relationship: key to accurate forecast is the choice of appropriate methods of data-mining and utilisation. In grey theory [2, 3], every stochastic process has grey variables changing at a certain amplitude and periodicity. An accumulated generating operation (AGO) is a basic method to render a grey process white. With accumulation, the changing trend becomes more apparent and the innate governing relationship is revealed. Based on this, a grey differential equation is built to describe the dynamic law of grey accumulation. Then, an inverse accumulated generating operation (IAGO) is used to reduce the whitening of this grey process. Of course, it is not merely the case that greater accumulation leads to better forecasts. Tien demonstrated that the k-AGO (k
>
1) data of the original series cannot be used as intermediate information in grey prediction model building [4]. In practical applications, the system generated by real data is considered grey: on that basis, grey models can be both effective past simulators and future predictors.In general, a grey model can be denoted by GM (m
,
n), where m is the order and n is the number of variables of the grey equation. GM (1, 1) is the most widely used grey prediction model with the simplest structure [2]. Its core principle is to regard real systems as generalised energy systems with exponentially changing trends under no external agency’s influence. Therefore, only when 1-AGO original data series are consistent with this exponential change, can accurate prediction by GM (1, 1) be achieved. Many scholars have made useful improvements to basic GM (1, 1) models [5–8]. In spite of this, it is not appropriate to use GM (1, 1) predictions based upon data mining from series which are themselves subject to external agents. The GM (1
,
n) model with n
-
1 relative factors acting as the associated series is a multiple, grey, prediction model [2, 3]. This model can fully exploit information provided the associated series unlike a GM (
1
,
1
) model which only includes information about the predicted series in its modelling process. So, from the point of view of provision of supplementary information, the prediction accuracy of GM (
1
,
n
), should, under these conditions, be higher than that of GM (1, 1). However, the solution of the whitening differential equation of GM (
1
,
n
) is inexact and on occasion wrong thereby producing significant practical forecasting errors [4]. The model therefore has found few real-world applications [9–11] to date.A grey prediction model with convolution integral GMC(
1
,
n
) proposed by Tien [12] is a new model to improve upon traditional GM (
1
,
n
) versions. The values modelled by GMC (1
,
n) are theoretically the exact solution of the traditional GM (
1
,
n
) model, and the grey control parameter like that of GM (1, 1) is introduced into the model. GMC (
1
,
n
) reduces to GM (1, 1) for n
=
1. The GMC (
1
,
n
) model improved the prediction accuracy of multiple grey models and has successfully been applied to different areas [12–14]. Tien proposed three improved models over the basic GMC (
1
,
n
) model to meet various application requirements. The three models are deterministic GMC (
1
,
n
) (DGDMC (1
,
n)) [15]; interval GDMC (
1
,
n
) (IGDMC (
1
,
n
)) [16]; and first pair-of-data GMC (
1
,
n
) (FGMC (
1
,
n
)) [17]. The first derivative of the 1-AGO data of each associated series is introduced into the DGDMC (1
,
n) model [15] to strengthen the indicative significance while evaluating its 1-AGO predicted data by a convolution integral. The IGDMC (
1
,
n
) model [16] draws lessons from the prediction method of its linear regression interval, extending the prediction method of DGDMC (1
,
n) from point to interval prediction. Certain components of the model may be removed to better satisfy hypothesis tests based upon system parameters. The building of the improved model FGMC (
1
,
n
) [17] only needs n
+
2 pairs of historical data: for the first pair of entries, the message was shown to be independent from the modelling results. Therefore, FGMC (1
,
n) is usually more satisfactory and stable than GMC (
1
,
n
) because it extracts the first pair of entries and bases its predictions on the message enshrined in the original series.Nowadays, GMC(
1
,
n
) and variants thereon are linear: the high-precision predictions are predicated upon the similarity in trend between the predicted variables and their influencing factors. At present, most applications meet this premise. For improved predictions, Wu and Chen [13] used grey relational analysis to analyse the similarities between predicted, and associated, series curves of differing periodicities. They then built GMC (
1
,
n
) using data from the largest period of grey relational analysis and obtained satisfactory predictions. However, the nonlinear interaction between the internal and external system factors is ubiquitous in reality: the linear structure of these GMC (
1
,
n
) models is not conducive to high-precision prediction of dynamic systems. Therefore, nonlinear processing to the forms of the relative factors on the right-hand side of the GMC (
1
,
n
) equation was applied here. A power exponent was added to the 1-AGO series of every associated series to reflect the nonlinear interaction of the associated series with that predicted. Although these newly added parameters were initially unknown, they were estimated and subjected to refinement by a least-squares method. After obtaining the time response function of the predicted series, the optimum values of parameters reflecting their nonlinear interactions were derived with regard to the aforementioned minimum error criterion. This improved GMC (
1
,
n
) model, a nonlinear grey prediction model with the addition of a convolution integral, is hereafter designated NGMC (
1
,
n
).In recent years, industrial SO2 emissions have become a major source of harmful air pollution in China as their industrial growth continues. To cater for these increased emissions, an enhanced understanding of likely future industrial SO2 emissions is needed. Though China’s industrial output values can be used as a basis to interpret growing industrial emissions, there is no a simultaneous increase in industrial emissions and industrial development [18]. Data also show that China’s SO2 emissions and industrial output values are not governed by a simple linear relationship but an uncertain nonlinear one. We present GMC (
1
,
n
) and NGMC (
1
,
n
) models where the influencing, causative, factors are industrial output; predicted series is China’s SO2 emissions. Through comparison of their modelling and prediction accuracy, the paper demonstrated the effectiveness of the proposed NGMC (
1
,
n
).The remainder of this paper is organized as follows. An introduction to the existing GMC(
1
,
n
) model and the modelling method of nonlinear grey model with convolution integral NGMC (1
,
n) proposed in this paper are given in Section 2. Section 3 demonstrates the effectiveness of the proposed NGMC (
1
,
n
) model by forecasting China’s SO2 emissions. Finally, the paper concludes with some comments in Section 4.


## 2. Modelling Method

### 2.1. The Existing GMC(
1
,
n
) Model

Assume that pairs of observationsX
(
0
)
=
(
X
1
(
0
)
,
X
2
(
0
)
,
…
,
X
n
(
0
)
) are available at time intervals of n
-
1 inputs X
2
(
0
)
,
X
3
(
0
)
,
…
,
X
n
(
0
) and an outputX
1
(
0
)from some dynamic system. The existing GMC (1
,
n) modelling process [12] is carried out as follows.Consider the following original predicted series:(1)
X
1
(
0
)
=
{
X
1
(
0
)
(
r
p
+
1
)
,
X
1
(
0
)
(
r
p
+
2
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
}
;
and the original associated series:
(2)
X
i
(
0
)
=
{
X
i
(
0
)
(
1
)
,
X
i
(
0
)
(
2
)
,
…
,
X
i
(
0
)
(
r
)
}
,
i
=
2,3
,
…
,
n
then the first order accumulated generation (1-AGO) data for X
1
(
0
)
,
X
2
(
0
)
,
…
,
X
n
(
0
) are given by (3) and (4), respectively:
(3)
X
1
(
1
)
(
r
p
+
k
)
=
∑
j
=
1
k
X
1
(
0
)
(
r
p
+
j
)
,
(4)
X
i
(
1
)
(
k
)
=
∑
j
=
1
k
X
i
(
0
)
(
j
)
,
i
=
2,3
,
…
,
n
.
The grey prediction model based on the predicted 1-AGO series:
(5)
X
1
(
1
)
=
{
X
1
(
1
)
(
r
p
+
1
)
,
X
1
(
1
)
(
r
p
+
2
)
,
…
,
X
1
(
1
)
(
r
p
+
r
)
}
and the associated 1-AGO series:
(6)
X
i
(
1
)
=
{
X
i
(
1
)
(
1
)
,
X
i
(
1
)
(
2
)
,
…
,
X
i
(
1
)
(
r
)
}
,
i
=
2,3
,
…
,
n
is given by the differential equation:
(7)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
+
b
1
X
1
(
1
)
(
r
p
+
t
)
=
b
2
X
2
(
1
)
(
t
)
+
b
3
X
3
(
1
)
(
t
)
+
⋯
+
b
n
X
n
(
1
)
(
t
)
+
u
,
1111111111111111111111111
.
t
=
1,2
,
…
,
r
f
,
where b
1
,
b
2
,
…
,
b
n,  u are parameters to be estimated and r is the data number used in model building; r
p is a delay period, and r
f is the number of entries to be forecast. Equation (7) is called the n-factor grey prediction model with convolution integral and is denoted by GMC (1
,
n) [12]; the 1 represents the first-order derivative of the 1-AGO series of X
1
(
1
); the n represents the total of n relative series introduced into the grey differential equation.The grey derivative for the first-order grey differential equation with 1-AGO is represented as follows:(8)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
lim
⁡
Δ
t
→
0
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
Δ
t
)
Δ
t
,
and when Δ
t
→
1
(9)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
Δ
X
1
(
1
)
(
r
p
+
t
)
Δ
t
=
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
1
)
=
X
1
(
0
)
(
r
p
+
t
)
.The background value of the grey derivatived
X
1
(
1
)
(
r
p
+
t
)
/
d
t is taken as the mean of X
1
(
1
)
(
r
p
+
t
) and X
1
(
1
)
(
r
p
+
t
-
1
), and those of the associated series X
i
(
1
)
(
t
) are also taken as the mean of X
i
(
1
)
(
t
) and X
i
(
1
)
(
t
-
1
) for i
=
2,3
,
…
,
n, respectively, in the determination of model parameters by GMC (1
,
n).The least-squares solution to the model parameters of GMC (1
,
n) [12] in (7) by t from 1 to r is
(10)
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
R
,
where (11)
B
=
[
-
X
1
(
1
)
(
r
p
+
1
)
+
X
1
(
1
)
(
r
p
+
2
)
2
X
2
(
1
)
(
1
)
+
X
2
(
1
)
(
2
)
2
X
3
(
1
)
(
1
)
+
X
3
(
1
)
(
2
)
2
⋯
X
n
(
1
)
(
1
)
+
X
n
(
1
)
(
2
)
2
1
-
X
1
(
1
)
(
r
p
+
2
)
+
X
1
(
1
)
(
r
p
+
3
)
2
X
2
(
1
)
(
2
)
+
X
2
(
1
)
(
3
)
2
X
3
(
1
)
(
2
)
+
X
3
(
1
)
(
3
)
2
⋯
X
n
(
1
)
(
2
)
+
X
n
(
1
)
(
3
)
2
1
⋯
⋯
⋯
⋯
⋯
⋯
-
X
1
(
1
)
(
r
p
+
r
-
1
)
+
X
1
(
1
)
(
r
p
+
r
)
2
X
2
(
1
)
(
r
-
1
)
+
X
2
(
1
)
(
r
)
2
X
3
(
1
)
(
r
-
1
)
+
X
3
(
1
)
(
r
)
2
⋯
X
n
(
1
)
(
r
-
1
)
+
X
n
(
1
)
(
r
)
2
1
]
,
(12)
Y
R
=
(
X
1
(
0
)
(
r
p
+
2
)
,
X
1
(
0
)
(
r
p
+
3
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
)
T
.In summary, the right-hand side of (7), the discrete function f
(
t
) [12] can be obtained as
(13)
f
(
t
)
=
b
2
X
2
(
1
)
(
t
)
+
b
3
X
3
(
1
)
(
t
)
+
⋯
+
b
n
X
n
(
1
)
(
t
)
+
u
,
1111111111111111111111111
t
=
1,2
,
…
,
r
f
.
The 1-AGO modelling values of the predicted series [12] can be derived with the initial condition X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
) as
(14)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
1
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
∫
1
t
e
-
b
1
(
t
-
τ
)
f
(
τ
)
d
τ
111111111111111111111111
t
=
1,2
,
…
,
r
f
.
The modelling values, X
^
1
(
1
)
(
r
p
+
t
),  t
=
1,2
,
…
,
r
f in (14) of the 1-AGO data of the predicted series [12] can be evaluated approximately by
(15)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(16)
X
^
1
(
1
)
(
r
p
+
t
)
≅
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
(
t
)
+
f
(
t
-
1
)
]
}
,
where u
(
t
-
2
) is the unit step function.Applying 1-IAGO to (16) yields the following modelled values together with the forecasts:
(17)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(18)
X
^
1
(
0
)
(
r
p
+
t
)
=
X
^
1
(
1
)
(
r
p
+
t
)
-
X
^
1
(
1
)
(
r
p
+
t
-
1
)
,
1111111111111111
t
=
2,3
,
…
,
r
f
.Assume the system parameters in (7) to be constants in the postsampling period and then by using the postsampling data, combined with the given data for the corresponding associated series, as a new input series, the corresponding forecasts or values of indirect measurement for the predicted series can be derived.


### 2.2. Nonlinear Grey Modelling with Convolution Integral: NGMC(
1
,
n
)

To adapt GMC (1
,
n) to reality’s ubiquitous nonlinearity, the improved model—NGMC (1
,
n)—with its convolution integral was developed and used.


#### 2.2.1. Representation of the New Model: NGMC(
1
,
n
)

Suppose that the grey prediction model based on the predicted 1-AGO series:(19)
X
1
(
1
)
=
{
X
1
(
1
)
(
r
p
+
1
)
,
X
1
(
1
)
(
r
p
+
2
)
,
…
,
X
1
(
1
)
(
r
p
+
r
)
}
and the associated 1-AGO series:
(20)
X
i
(
1
)
=
{
X
i
(
1
)
(
1
)
,
X
i
(
1
)
(
2
)
,
…
,
X
i
(
1
)
(
r
)
}
,
i
=
2,3
,
…
,
n
is given by the differential equation:
(21)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
+
b
1
X
1
(
1
)
(
r
p
+
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
111111111111111111111111111.111111111
t
=
1,2
,
…
,
r
f
,
where b
1
,
b
2
,
…
,
b
n, u, β
2
,
β
3
,
…
,
β
n are parameters to be estimated, r is the data number used in model building, r
p is a delay period, and r
f the number of entries to be forecast. Equation (21) denotes the n-factor nonlinear grey prediction model with the addition of a convolution integral (NGMC (1
,
n)), where 1 represents the first-order derivative of the 1-AGO series of X
1
(
1
) and n represents the total numberrelative series introduced into the grey differential equation.Compared with GMC(
1
,
n
),  n
-
1unknown parameters β
2
,
β
3
,
…
,
β
n are introduced into this NGMC (1
,
n) model. They are taken as the n
-
1 power exponents of the predicted variables’ relevant factors to reflect the effect of these upon nonlinear system behaviours and interactions.Whenβ
2
=
β
3
=
⋯
=
β
n
=
1, (20) is reduced to a traditional GMC (
1
,
n
) model [12].
Whenβ
2
=
β
3
=
⋯
=
β
n
=
1 and u
=
0, (20) reduces to a traditional GM (1
,
n) model [2, 4].
Whenn
=
1 or β
2
=
β
3
=
⋯
=
β
n
=
0, (20) reduces to a traditional GM (1, 1) model [2].


#### 2.2.2. The Evaluation of System Parametersb
1
,
b
2
,
…
,
b
n, u, and f
N
(
t
)

To evaluate all parameters in the model, the power exponentsβ
2
,
β
3
,
…
,
β
n of X
2
(
1
)
,
X
3
(
1
)
,
…
,
X
n
(
1
) are assumed to be known; by applying a method similar to that of GMC (1
,
n), estimates of the corresponding parameters b
1
,
b
2
,
…
,
b
n, u are made. After obtaining the time response function of X
^
1
(
0
), the corresponding optimisation algorithm is applied to solve for β
2
,
β
3
,
…
,
β
n.The grey derivative for the first-order grey differential equation with 1-AGO is conventionally represented by(22)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
1
)
=
X
1
(
0
)
(
r
p
+
t
)
.
The background value of the grey derivative d
X
1
(
1
)
(
r
p
+
t
)
/
d
t is taken as the means of X
1
(
1
)
(
r
p
+
t
) and X
1
(
1
)
(
r
p
+
t
-
1
), and those of the associated series X
i
(
1
)
(
t
) are also taken as the means of X
i
(
1
)
(
t
) and X
i
(
1
)
(
t
-
1
) for i
=
2,3
,
…
,
n, respectively, when determining the model parameters by NGMC (
1
,
n
).The least-squares solution to the model parameters of NGMC (1
,
n) in (20) by t from 1 to r is
(23)
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
R
,
where(24)
B
=
[
-
X
1
(
1
)
(
r
p
+
1
)
+
X
1
(
1
)
(
r
p
+
2
)
2
[
X
2
(
1
)
(
1
)
+
X
2
(
1
)
(
2
)
2
]
β
2
[
X
3
(
1
)
(
1
)
+
X
3
(
1
)
(
2
)
2
]
β
3
⋯
[
X
n
(
1
)
(
1
)
+
X
n
(
1
)
(
2
)
2
]
β
n
1
-
X
1
(
1
)
(
r
p
+
2
)
+
X
1
(
1
)
(
r
p
+
3
)
2
[
X
2
(
1
)
(
2
)
+
X
2
(
1
)
(
3
)
2
]
β
2
[
X
3
(
1
)
(
2
)
+
X
3
(
1
)
(
3
)
2
]
β
3
⋯
[
X
n
(
1
)
(
2
)
+
X
n
(
1
)
(
3
)
2
]
β
n
1
⋯
⋯
⋯
⋯
⋯
⋯
-
X
1
(
1
)
(
r
p
+
r
-
1
)
+
X
1
(
1
)
(
r
p
+
r
)
2
[
X
2
(
1
)
(
r
-
1
)
+
X
2
(
1
)
(
r
)
2
]
β
2
[
X
3
(
1
)
(
r
-
1
)
+
X
3
(
1
)
(
r
)
2
]
β
3
⋯
[
X
n
(
1
)
(
r
-
1
)
+
X
n
(
1
)
(
r
)
2
]
β
n
1
]
,
(25)
Y
R
=
(
X
1
(
0
)
(
r
p
+
2
)
,
X
1
(
0
)
(
r
p
+
3
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
)
T
.In summary, the right-hand side of (21), the discrete function f
(
t
), can be obtained as
(26)
f
N
(
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
111111111
t
=
1,2
,
…
,
r
f
.As seen in (26), f
N
(
t
) is a nonlinear function on X
i
(
1
),  i
=
2,3
,
…
,
n, but it, in essence, is still linear. If the new variable Y
i
(
1
)
(
t
)
=
[
X
i
(
1
)
(
t
)
]
β
i,  i
=
2,3
,
…
,
n, f
N
(
t
) can be written as
(27)
f
N
(
t
)
=
b
2
Y
2
(
1
)
(
t
)
+
b
3
Y
3
(
1
)
(
t
)
+
⋯
+
b
n
Y
n
(
1
)
(
t
)
+
u
,
t
=
1,2
,
…
,
r
f
;
therefore, (27) remains a linear function, essentially the same as (13).


#### 2.2.3. The Determination of Unit Impulse Response Functionh
(
t
)

The unit impulse response functionh
(
t
) of the system characterised by (21) can be derived by Laplace transform. From (21):
(28)
d
X
1
(
1
)
(
t
)
d
t
+
b
1
X
1
(
1
)
(
t
)
=
δ
(
t
)
,
where δ
(
t
) is the unit impulse function to determine the corresponding unit impulse response function. Applying a Laplace transform to (28) with the initial condition X
1
(
1
)
(
1
)
=
0 gives
(29)
s
X
-
1
(
1
)
(
s
)
+
b
1
X
-
1
(
1
)
(
s
)
=
1
or
(30)
X
-
1
(
1
)
(
s
)
=
1
s
+
b
1
.
The inverse transform of X
-
1
(
1
)
(
s
) is
(31)
X
1
(
1
)
(
t
)
=
e
-
b
1
t
.
That is, the unit impulse response function h
(
t
) of the system is
(32)
h
(
t
)
=
e
-
b
1
t
.
The process continues by outlining the evaluation of X
^
1
(
0
).


#### 2.2.4. The Evaluation ofX
^
1
(
0
)

The 1-AGO modelling values of the predicted series can be derived with the initial conditionX
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
) as
(33)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
1
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
∫
1
t
e
-
b
1
(
t
-
τ
)
f
N
(
τ
)
d
τ
,
111111
t
=
1,2
,
…
,
r
f
.
The second term on the right-hand side of (33) can be evaluated approximately by two-point Gaussian numerical integration. Thus, the modelling values in (33) of the 1-AGO data of the predicted series can be evaluated approximately by
(34)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(35)
X
^
1
(
1
)
(
r
p
+
t
)
≅
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
N
(
t
)
+
f
N
(
t
-
1
)
]
}
,
where u
(
t
-
2
) is the unit step function.Applying 1-IAGO to (35) produced the following modelling values and their forecasts:
(36)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(37)
X
^
1
(
0
)
(
r
p
+
t
)
=
X
^
1
(
1
)
(
r
p
+
t
)
-
X
^
1
(
1
)
(
r
p
+
t
-
1
)
,
11111111111111111
t
=
2,3
,
…
,
r
f
.


#### 2.2.5. An Overall Measure of Accuracy for Forecasts

To evaluate forecast performance, Tien’s standard [12] test using the root mean squared percentage error (RMSPE) for the priori-sample period (RMSPEPR) and postsample periods (RMSPEPO), respectively, is used. Generally, the RMSPEPR and RMSPEPO are defined as
(38)
RMSPEPR
=
1
r
∑
t
=
r
p
+
1
r
p
+
r
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
,
(39)
RMSPEPO
=
1
r
f
∑
t
=
r
p
+
r
+
1
r
p
+
r
+
r
f
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
.


#### 2.2.6. The Determination of Power Exponents

The discussion above is proposed under the condition thatβ
i,  i
=
2,3
,
…
,
n is known, whereas, in reality, we cannot really know the values of them that are unknown. This work suggests that the minimisation of the root mean squared percentage error for the presample period should be the objective, and the unknown parameters β
i
=
1,  i
=
2,3
,
…
,
n are solved by building the optimisation model below.Consider the following:(40)
Min
⁡
β
i
RMSPEPR
=
1
r
∑
t
=
r
p
+
1
r
p
+
r
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
,
i
=
2,3
,
…
,
n
s
.
t
.
{
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
N
(
t
)
+
f
N
(
t
-
1
)
]
}
,
t
=
2,3
,
…
,
r
f
N
(
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
t
=
1,2
,
…
,
r
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
RThe optimisation problem above can be solved by proprietary software. Onceβ
i
=
1, i
=
2,3
,
…
,
n is confirmed, the structure parameters b
1
,
b
2
,
…
,
b
n, u of the model are also determined. Then, the discrete function f
N
(
t
) can be derived by back-substituting the derived results into (26): final predictions can then be obtained by applying f
N
(
t
) to (35) to (37).


## 2.1. The Existing GMC(
1
,
n
) Model

Assume that pairs of observationsX
(
0
)
=
(
X
1
(
0
)
,
X
2
(
0
)
,
…
,
X
n
(
0
)
) are available at time intervals of n
-
1 inputs X
2
(
0
)
,
X
3
(
0
)
,
…
,
X
n
(
0
) and an outputX
1
(
0
)from some dynamic system. The existing GMC (1
,
n) modelling process [12] is carried out as follows.Consider the following original predicted series:(1)
X
1
(
0
)
=
{
X
1
(
0
)
(
r
p
+
1
)
,
X
1
(
0
)
(
r
p
+
2
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
}
;
and the original associated series:
(2)
X
i
(
0
)
=
{
X
i
(
0
)
(
1
)
,
X
i
(
0
)
(
2
)
,
…
,
X
i
(
0
)
(
r
)
}
,
i
=
2,3
,
…
,
n
then the first order accumulated generation (1-AGO) data for X
1
(
0
)
,
X
2
(
0
)
,
…
,
X
n
(
0
) are given by (3) and (4), respectively:
(3)
X
1
(
1
)
(
r
p
+
k
)
=
∑
j
=
1
k
X
1
(
0
)
(
r
p
+
j
)
,
(4)
X
i
(
1
)
(
k
)
=
∑
j
=
1
k
X
i
(
0
)
(
j
)
,
i
=
2,3
,
…
,
n
.
The grey prediction model based on the predicted 1-AGO series:
(5)
X
1
(
1
)
=
{
X
1
(
1
)
(
r
p
+
1
)
,
X
1
(
1
)
(
r
p
+
2
)
,
…
,
X
1
(
1
)
(
r
p
+
r
)
}
and the associated 1-AGO series:
(6)
X
i
(
1
)
=
{
X
i
(
1
)
(
1
)
,
X
i
(
1
)
(
2
)
,
…
,
X
i
(
1
)
(
r
)
}
,
i
=
2,3
,
…
,
n
is given by the differential equation:
(7)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
+
b
1
X
1
(
1
)
(
r
p
+
t
)
=
b
2
X
2
(
1
)
(
t
)
+
b
3
X
3
(
1
)
(
t
)
+
⋯
+
b
n
X
n
(
1
)
(
t
)
+
u
,
1111111111111111111111111
.
t
=
1,2
,
…
,
r
f
,
where b
1
,
b
2
,
…
,
b
n,  u are parameters to be estimated and r is the data number used in model building; r
p is a delay period, and r
f is the number of entries to be forecast. Equation (7) is called the n-factor grey prediction model with convolution integral and is denoted by GMC (1
,
n) [12]; the 1 represents the first-order derivative of the 1-AGO series of X
1
(
1
); the n represents the total of n relative series introduced into the grey differential equation.The grey derivative for the first-order grey differential equation with 1-AGO is represented as follows:(8)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
lim
⁡
Δ
t
→
0
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
Δ
t
)
Δ
t
,
and when Δ
t
→
1
(9)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
Δ
X
1
(
1
)
(
r
p
+
t
)
Δ
t
=
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
1
)
=
X
1
(
0
)
(
r
p
+
t
)
.The background value of the grey derivatived
X
1
(
1
)
(
r
p
+
t
)
/
d
t is taken as the mean of X
1
(
1
)
(
r
p
+
t
) and X
1
(
1
)
(
r
p
+
t
-
1
), and those of the associated series X
i
(
1
)
(
t
) are also taken as the mean of X
i
(
1
)
(
t
) and X
i
(
1
)
(
t
-
1
) for i
=
2,3
,
…
,
n, respectively, in the determination of model parameters by GMC (1
,
n).The least-squares solution to the model parameters of GMC (1
,
n) [12] in (7) by t from 1 to r is
(10)
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
R
,
where (11)
B
=
[
-
X
1
(
1
)
(
r
p
+
1
)
+
X
1
(
1
)
(
r
p
+
2
)
2
X
2
(
1
)
(
1
)
+
X
2
(
1
)
(
2
)
2
X
3
(
1
)
(
1
)
+
X
3
(
1
)
(
2
)
2
⋯
X
n
(
1
)
(
1
)
+
X
n
(
1
)
(
2
)
2
1
-
X
1
(
1
)
(
r
p
+
2
)
+
X
1
(
1
)
(
r
p
+
3
)
2
X
2
(
1
)
(
2
)
+
X
2
(
1
)
(
3
)
2
X
3
(
1
)
(
2
)
+
X
3
(
1
)
(
3
)
2
⋯
X
n
(
1
)
(
2
)
+
X
n
(
1
)
(
3
)
2
1
⋯
⋯
⋯
⋯
⋯
⋯
-
X
1
(
1
)
(
r
p
+
r
-
1
)
+
X
1
(
1
)
(
r
p
+
r
)
2
X
2
(
1
)
(
r
-
1
)
+
X
2
(
1
)
(
r
)
2
X
3
(
1
)
(
r
-
1
)
+
X
3
(
1
)
(
r
)
2
⋯
X
n
(
1
)
(
r
-
1
)
+
X
n
(
1
)
(
r
)
2
1
]
,
(12)
Y
R
=
(
X
1
(
0
)
(
r
p
+
2
)
,
X
1
(
0
)
(
r
p
+
3
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
)
T
.In summary, the right-hand side of (7), the discrete function f
(
t
) [12] can be obtained as
(13)
f
(
t
)
=
b
2
X
2
(
1
)
(
t
)
+
b
3
X
3
(
1
)
(
t
)
+
⋯
+
b
n
X
n
(
1
)
(
t
)
+
u
,
1111111111111111111111111
t
=
1,2
,
…
,
r
f
.
The 1-AGO modelling values of the predicted series [12] can be derived with the initial condition X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
) as
(14)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
1
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
∫
1
t
e
-
b
1
(
t
-
τ
)
f
(
τ
)
d
τ
111111111111111111111111
t
=
1,2
,
…
,
r
f
.
The modelling values, X
^
1
(
1
)
(
r
p
+
t
),  t
=
1,2
,
…
,
r
f in (14) of the 1-AGO data of the predicted series [12] can be evaluated approximately by
(15)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(16)
X
^
1
(
1
)
(
r
p
+
t
)
≅
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
(
t
)
+
f
(
t
-
1
)
]
}
,
where u
(
t
-
2
) is the unit step function.Applying 1-IAGO to (16) yields the following modelled values together with the forecasts:
(17)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(18)
X
^
1
(
0
)
(
r
p
+
t
)
=
X
^
1
(
1
)
(
r
p
+
t
)
-
X
^
1
(
1
)
(
r
p
+
t
-
1
)
,
1111111111111111
t
=
2,3
,
…
,
r
f
.Assume the system parameters in (7) to be constants in the postsampling period and then by using the postsampling data, combined with the given data for the corresponding associated series, as a new input series, the corresponding forecasts or values of indirect measurement for the predicted series can be derived.


## 2.2. Nonlinear Grey Modelling with Convolution Integral: NGMC(
1
,
n
)

To adapt GMC (1
,
n) to reality’s ubiquitous nonlinearity, the improved model—NGMC (1
,
n)—with its convolution integral was developed and used.


### 2.2.1. Representation of the New Model: NGMC(
1
,
n
)

Suppose that the grey prediction model based on the predicted 1-AGO series:(19)
X
1
(
1
)
=
{
X
1
(
1
)
(
r
p
+
1
)
,
X
1
(
1
)
(
r
p
+
2
)
,
…
,
X
1
(
1
)
(
r
p
+
r
)
}
and the associated 1-AGO series:
(20)
X
i
(
1
)
=
{
X
i
(
1
)
(
1
)
,
X
i
(
1
)
(
2
)
,
…
,
X
i
(
1
)
(
r
)
}
,
i
=
2,3
,
…
,
n
is given by the differential equation:
(21)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
+
b
1
X
1
(
1
)
(
r
p
+
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
111111111111111111111111111.111111111
t
=
1,2
,
…
,
r
f
,
where b
1
,
b
2
,
…
,
b
n, u, β
2
,
β
3
,
…
,
β
n are parameters to be estimated, r is the data number used in model building, r
p is a delay period, and r
f the number of entries to be forecast. Equation (21) denotes the n-factor nonlinear grey prediction model with the addition of a convolution integral (NGMC (1
,
n)), where 1 represents the first-order derivative of the 1-AGO series of X
1
(
1
) and n represents the total numberrelative series introduced into the grey differential equation.Compared with GMC(
1
,
n
),  n
-
1unknown parameters β
2
,
β
3
,
…
,
β
n are introduced into this NGMC (1
,
n) model. They are taken as the n
-
1 power exponents of the predicted variables’ relevant factors to reflect the effect of these upon nonlinear system behaviours and interactions.Whenβ
2
=
β
3
=
⋯
=
β
n
=
1, (20) is reduced to a traditional GMC (
1
,
n
) model [12].
Whenβ
2
=
β
3
=
⋯
=
β
n
=
1 and u
=
0, (20) reduces to a traditional GM (1
,
n) model [2, 4].
Whenn
=
1 or β
2
=
β
3
=
⋯
=
β
n
=
0, (20) reduces to a traditional GM (1, 1) model [2].


### 2.2.2. The Evaluation of System Parametersb
1
,
b
2
,
…
,
b
n, u, and f
N
(
t
)

To evaluate all parameters in the model, the power exponentsβ
2
,
β
3
,
…
,
β
n of X
2
(
1
)
,
X
3
(
1
)
,
…
,
X
n
(
1
) are assumed to be known; by applying a method similar to that of GMC (1
,
n), estimates of the corresponding parameters b
1
,
b
2
,
…
,
b
n, u are made. After obtaining the time response function of X
^
1
(
0
), the corresponding optimisation algorithm is applied to solve for β
2
,
β
3
,
…
,
β
n.The grey derivative for the first-order grey differential equation with 1-AGO is conventionally represented by(22)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
1
)
=
X
1
(
0
)
(
r
p
+
t
)
.
The background value of the grey derivative d
X
1
(
1
)
(
r
p
+
t
)
/
d
t is taken as the means of X
1
(
1
)
(
r
p
+
t
) and X
1
(
1
)
(
r
p
+
t
-
1
), and those of the associated series X
i
(
1
)
(
t
) are also taken as the means of X
i
(
1
)
(
t
) and X
i
(
1
)
(
t
-
1
) for i
=
2,3
,
…
,
n, respectively, when determining the model parameters by NGMC (
1
,
n
).The least-squares solution to the model parameters of NGMC (1
,
n) in (20) by t from 1 to r is
(23)
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
R
,
where(24)
B
=
[
-
X
1
(
1
)
(
r
p
+
1
)
+
X
1
(
1
)
(
r
p
+
2
)
2
[
X
2
(
1
)
(
1
)
+
X
2
(
1
)
(
2
)
2
]
β
2
[
X
3
(
1
)
(
1
)
+
X
3
(
1
)
(
2
)
2
]
β
3
⋯
[
X
n
(
1
)
(
1
)
+
X
n
(
1
)
(
2
)
2
]
β
n
1
-
X
1
(
1
)
(
r
p
+
2
)
+
X
1
(
1
)
(
r
p
+
3
)
2
[
X
2
(
1
)
(
2
)
+
X
2
(
1
)
(
3
)
2
]
β
2
[
X
3
(
1
)
(
2
)
+
X
3
(
1
)
(
3
)
2
]
β
3
⋯
[
X
n
(
1
)
(
2
)
+
X
n
(
1
)
(
3
)
2
]
β
n
1
⋯
⋯
⋯
⋯
⋯
⋯
-
X
1
(
1
)
(
r
p
+
r
-
1
)
+
X
1
(
1
)
(
r
p
+
r
)
2
[
X
2
(
1
)
(
r
-
1
)
+
X
2
(
1
)
(
r
)
2
]
β
2
[
X
3
(
1
)
(
r
-
1
)
+
X
3
(
1
)
(
r
)
2
]
β
3
⋯
[
X
n
(
1
)
(
r
-
1
)
+
X
n
(
1
)
(
r
)
2
]
β
n
1
]
,
(25)
Y
R
=
(
X
1
(
0
)
(
r
p
+
2
)
,
X
1
(
0
)
(
r
p
+
3
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
)
T
.In summary, the right-hand side of (21), the discrete function f
(
t
), can be obtained as
(26)
f
N
(
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
111111111
t
=
1,2
,
…
,
r
f
.As seen in (26), f
N
(
t
) is a nonlinear function on X
i
(
1
),  i
=
2,3
,
…
,
n, but it, in essence, is still linear. If the new variable Y
i
(
1
)
(
t
)
=
[
X
i
(
1
)
(
t
)
]
β
i,  i
=
2,3
,
…
,
n, f
N
(
t
) can be written as
(27)
f
N
(
t
)
=
b
2
Y
2
(
1
)
(
t
)
+
b
3
Y
3
(
1
)
(
t
)
+
⋯
+
b
n
Y
n
(
1
)
(
t
)
+
u
,
t
=
1,2
,
…
,
r
f
;
therefore, (27) remains a linear function, essentially the same as (13).


### 2.2.3. The Determination of Unit Impulse Response Functionh
(
t
)

The unit impulse response functionh
(
t
) of the system characterised by (21) can be derived by Laplace transform. From (21):
(28)
d
X
1
(
1
)
(
t
)
d
t
+
b
1
X
1
(
1
)
(
t
)
=
δ
(
t
)
,
where δ
(
t
) is the unit impulse function to determine the corresponding unit impulse response function. Applying a Laplace transform to (28) with the initial condition X
1
(
1
)
(
1
)
=
0 gives
(29)
s
X
-
1
(
1
)
(
s
)
+
b
1
X
-
1
(
1
)
(
s
)
=
1
or
(30)
X
-
1
(
1
)
(
s
)
=
1
s
+
b
1
.
The inverse transform of X
-
1
(
1
)
(
s
) is
(31)
X
1
(
1
)
(
t
)
=
e
-
b
1
t
.
That is, the unit impulse response function h
(
t
) of the system is
(32)
h
(
t
)
=
e
-
b
1
t
.
The process continues by outlining the evaluation of X
^
1
(
0
).


### 2.2.4. The Evaluation ofX
^
1
(
0
)

The 1-AGO modelling values of the predicted series can be derived with the initial conditionX
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
) as
(33)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
1
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
∫
1
t
e
-
b
1
(
t
-
τ
)
f
N
(
τ
)
d
τ
,
111111
t
=
1,2
,
…
,
r
f
.
The second term on the right-hand side of (33) can be evaluated approximately by two-point Gaussian numerical integration. Thus, the modelling values in (33) of the 1-AGO data of the predicted series can be evaluated approximately by
(34)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(35)
X
^
1
(
1
)
(
r
p
+
t
)
≅
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
N
(
t
)
+
f
N
(
t
-
1
)
]
}
,
where u
(
t
-
2
) is the unit step function.Applying 1-IAGO to (35) produced the following modelling values and their forecasts:
(36)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(37)
X
^
1
(
0
)
(
r
p
+
t
)
=
X
^
1
(
1
)
(
r
p
+
t
)
-
X
^
1
(
1
)
(
r
p
+
t
-
1
)
,
11111111111111111
t
=
2,3
,
…
,
r
f
.


### 2.2.5. An Overall Measure of Accuracy for Forecasts

To evaluate forecast performance, Tien’s standard [12] test using the root mean squared percentage error (RMSPE) for the priori-sample period (RMSPEPR) and postsample periods (RMSPEPO), respectively, is used. Generally, the RMSPEPR and RMSPEPO are defined as
(38)
RMSPEPR
=
1
r
∑
t
=
r
p
+
1
r
p
+
r
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
,
(39)
RMSPEPO
=
1
r
f
∑
t
=
r
p
+
r
+
1
r
p
+
r
+
r
f
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
.


### 2.2.6. The Determination of Power Exponents

The discussion above is proposed under the condition thatβ
i,  i
=
2,3
,
…
,
n is known, whereas, in reality, we cannot really know the values of them that are unknown. This work suggests that the minimisation of the root mean squared percentage error for the presample period should be the objective, and the unknown parameters β
i
=
1,  i
=
2,3
,
…
,
n are solved by building the optimisation model below.Consider the following:(40)
Min
⁡
β
i
RMSPEPR
=
1
r
∑
t
=
r
p
+
1
r
p
+
r
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
,
i
=
2,3
,
…
,
n
s
.
t
.
{
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
N
(
t
)
+
f
N
(
t
-
1
)
]
}
,
t
=
2,3
,
…
,
r
f
N
(
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
t
=
1,2
,
…
,
r
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
RThe optimisation problem above can be solved by proprietary software. Onceβ
i
=
1, i
=
2,3
,
…
,
n is confirmed, the structure parameters b
1
,
b
2
,
…
,
b
n, u of the model are also determined. Then, the discrete function f
N
(
t
) can be derived by back-substituting the derived results into (26): final predictions can then be obtained by applying f
N
(
t
) to (35) to (37).


## 2.2.1. Representation of the New Model: NGMC(
1
,
n
)

Suppose that the grey prediction model based on the predicted 1-AGO series:(19)
X
1
(
1
)
=
{
X
1
(
1
)
(
r
p
+
1
)
,
X
1
(
1
)
(
r
p
+
2
)
,
…
,
X
1
(
1
)
(
r
p
+
r
)
}
and the associated 1-AGO series:
(20)
X
i
(
1
)
=
{
X
i
(
1
)
(
1
)
,
X
i
(
1
)
(
2
)
,
…
,
X
i
(
1
)
(
r
)
}
,
i
=
2,3
,
…
,
n
is given by the differential equation:
(21)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
+
b
1
X
1
(
1
)
(
r
p
+
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
111111111111111111111111111.111111111
t
=
1,2
,
…
,
r
f
,
where b
1
,
b
2
,
…
,
b
n, u, β
2
,
β
3
,
…
,
β
n are parameters to be estimated, r is the data number used in model building, r
p is a delay period, and r
f the number of entries to be forecast. Equation (21) denotes the n-factor nonlinear grey prediction model with the addition of a convolution integral (NGMC (1
,
n)), where 1 represents the first-order derivative of the 1-AGO series of X
1
(
1
) and n represents the total numberrelative series introduced into the grey differential equation.Compared with GMC(
1
,
n
),  n
-
1unknown parameters β
2
,
β
3
,
…
,
β
n are introduced into this NGMC (1
,
n) model. They are taken as the n
-
1 power exponents of the predicted variables’ relevant factors to reflect the effect of these upon nonlinear system behaviours and interactions.Whenβ
2
=
β
3
=
⋯
=
β
n
=
1, (20) is reduced to a traditional GMC (
1
,
n
) model [12].
Whenβ
2
=
β
3
=
⋯
=
β
n
=
1 and u
=
0, (20) reduces to a traditional GM (1
,
n) model [2, 4].
Whenn
=
1 or β
2
=
β
3
=
⋯
=
β
n
=
0, (20) reduces to a traditional GM (1, 1) model [2].


## 2.2.2. The Evaluation of System Parametersb
1
,
b
2
,
…
,
b
n, u, and f
N
(
t
)

To evaluate all parameters in the model, the power exponentsβ
2
,
β
3
,
…
,
β
n of X
2
(
1
)
,
X
3
(
1
)
,
…
,
X
n
(
1
) are assumed to be known; by applying a method similar to that of GMC (1
,
n), estimates of the corresponding parameters b
1
,
b
2
,
…
,
b
n, u are made. After obtaining the time response function of X
^
1
(
0
), the corresponding optimisation algorithm is applied to solve for β
2
,
β
3
,
…
,
β
n.The grey derivative for the first-order grey differential equation with 1-AGO is conventionally represented by(22)
d
X
1
(
1
)
(
r
p
+
t
)
d
t
=
X
1
(
1
)
(
r
p
+
t
)
-
X
1
(
1
)
(
r
p
+
t
-
1
)
=
X
1
(
0
)
(
r
p
+
t
)
.
The background value of the grey derivative d
X
1
(
1
)
(
r
p
+
t
)
/
d
t is taken as the means of X
1
(
1
)
(
r
p
+
t
) and X
1
(
1
)
(
r
p
+
t
-
1
), and those of the associated series X
i
(
1
)
(
t
) are also taken as the means of X
i
(
1
)
(
t
) and X
i
(
1
)
(
t
-
1
) for i
=
2,3
,
…
,
n, respectively, when determining the model parameters by NGMC (
1
,
n
).The least-squares solution to the model parameters of NGMC (1
,
n) in (20) by t from 1 to r is
(23)
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
R
,
where(24)
B
=
[
-
X
1
(
1
)
(
r
p
+
1
)
+
X
1
(
1
)
(
r
p
+
2
)
2
[
X
2
(
1
)
(
1
)
+
X
2
(
1
)
(
2
)
2
]
β
2
[
X
3
(
1
)
(
1
)
+
X
3
(
1
)
(
2
)
2
]
β
3
⋯
[
X
n
(
1
)
(
1
)
+
X
n
(
1
)
(
2
)
2
]
β
n
1
-
X
1
(
1
)
(
r
p
+
2
)
+
X
1
(
1
)
(
r
p
+
3
)
2
[
X
2
(
1
)
(
2
)
+
X
2
(
1
)
(
3
)
2
]
β
2
[
X
3
(
1
)
(
2
)
+
X
3
(
1
)
(
3
)
2
]
β
3
⋯
[
X
n
(
1
)
(
2
)
+
X
n
(
1
)
(
3
)
2
]
β
n
1
⋯
⋯
⋯
⋯
⋯
⋯
-
X
1
(
1
)
(
r
p
+
r
-
1
)
+
X
1
(
1
)
(
r
p
+
r
)
2
[
X
2
(
1
)
(
r
-
1
)
+
X
2
(
1
)
(
r
)
2
]
β
2
[
X
3
(
1
)
(
r
-
1
)
+
X
3
(
1
)
(
r
)
2
]
β
3
⋯
[
X
n
(
1
)
(
r
-
1
)
+
X
n
(
1
)
(
r
)
2
]
β
n
1
]
,
(25)
Y
R
=
(
X
1
(
0
)
(
r
p
+
2
)
,
X
1
(
0
)
(
r
p
+
3
)
,
…
,
X
1
(
0
)
(
r
p
+
r
)
)
T
.In summary, the right-hand side of (21), the discrete function f
(
t
), can be obtained as
(26)
f
N
(
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
111111111
t
=
1,2
,
…
,
r
f
.As seen in (26), f
N
(
t
) is a nonlinear function on X
i
(
1
),  i
=
2,3
,
…
,
n, but it, in essence, is still linear. If the new variable Y
i
(
1
)
(
t
)
=
[
X
i
(
1
)
(
t
)
]
β
i,  i
=
2,3
,
…
,
n, f
N
(
t
) can be written as
(27)
f
N
(
t
)
=
b
2
Y
2
(
1
)
(
t
)
+
b
3
Y
3
(
1
)
(
t
)
+
⋯
+
b
n
Y
n
(
1
)
(
t
)
+
u
,
t
=
1,2
,
…
,
r
f
;
therefore, (27) remains a linear function, essentially the same as (13).


## 2.2.3. The Determination of Unit Impulse Response Functionh
(
t
)

The unit impulse response functionh
(
t
) of the system characterised by (21) can be derived by Laplace transform. From (21):
(28)
d
X
1
(
1
)
(
t
)
d
t
+
b
1
X
1
(
1
)
(
t
)
=
δ
(
t
)
,
where δ
(
t
) is the unit impulse function to determine the corresponding unit impulse response function. Applying a Laplace transform to (28) with the initial condition X
1
(
1
)
(
1
)
=
0 gives
(29)
s
X
-
1
(
1
)
(
s
)
+
b
1
X
-
1
(
1
)
(
s
)
=
1
or
(30)
X
-
1
(
1
)
(
s
)
=
1
s
+
b
1
.
The inverse transform of X
-
1
(
1
)
(
s
) is
(31)
X
1
(
1
)
(
t
)
=
e
-
b
1
t
.
That is, the unit impulse response function h
(
t
) of the system is
(32)
h
(
t
)
=
e
-
b
1
t
.
The process continues by outlining the evaluation of X
^
1
(
0
).


## 2.2.4. The Evaluation ofX
^
1
(
0
)

The 1-AGO modelling values of the predicted series can be derived with the initial conditionX
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
) as
(33)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
1
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
∫
1
t
e
-
b
1
(
t
-
τ
)
f
N
(
τ
)
d
τ
,
111111
t
=
1,2
,
…
,
r
f
.
The second term on the right-hand side of (33) can be evaluated approximately by two-point Gaussian numerical integration. Thus, the modelling values in (33) of the 1-AGO data of the predicted series can be evaluated approximately by
(34)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(35)
X
^
1
(
1
)
(
r
p
+
t
)
≅
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
N
(
t
)
+
f
N
(
t
-
1
)
]
}
,
where u
(
t
-
2
) is the unit step function.Applying 1-IAGO to (35) produced the following modelling values and their forecasts:
(36)
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
,
(37)
X
^
1
(
0
)
(
r
p
+
t
)
=
X
^
1
(
1
)
(
r
p
+
t
)
-
X
^
1
(
1
)
(
r
p
+
t
-
1
)
,
11111111111111111
t
=
2,3
,
…
,
r
f
.


## 2.2.5. An Overall Measure of Accuracy for Forecasts

To evaluate forecast performance, Tien’s standard [12] test using the root mean squared percentage error (RMSPE) for the priori-sample period (RMSPEPR) and postsample periods (RMSPEPO), respectively, is used. Generally, the RMSPEPR and RMSPEPO are defined as
(38)
RMSPEPR
=
1
r
∑
t
=
r
p
+
1
r
p
+
r
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
,
(39)
RMSPEPO
=
1
r
f
∑
t
=
r
p
+
r
+
1
r
p
+
r
+
r
f
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
.


## 2.2.6. The Determination of Power Exponents

The discussion above is proposed under the condition thatβ
i,  i
=
2,3
,
…
,
n is known, whereas, in reality, we cannot really know the values of them that are unknown. This work suggests that the minimisation of the root mean squared percentage error for the presample period should be the objective, and the unknown parameters β
i
=
1,  i
=
2,3
,
…
,
n are solved by building the optimisation model below.Consider the following:(40)
Min
⁡
β
i
RMSPEPR
=
1
r
∑
t
=
r
p
+
1
r
p
+
r
[
X
^
1
(
0
)
(
t
)
-
X
1
(
0
)
(
t
)
]
2
[
X
1
(
0
)
(
t
)
]
2
×
100
%
,
i
=
2,3
,
…
,
n
s
.
t
.
{
X
^
1
(
1
)
(
r
p
+
1
)
=
X
1
(
1
)
(
r
p
+
1
)
=
X
1
(
0
)
(
r
p
+
1
)
X
^
1
(
1
)
(
r
p
+
t
)
=
X
1
(
0
)
(
r
p
+
1
)
e
-
b
1
(
t
-
1
)
+
u
(
t
-
2
)
×
{
∑
k
=
2
t
e
-
b
1
(
t
-
k
+
(
1
/
2
)
)
·
1
2
[
f
N
(
t
)
+
f
N
(
t
-
1
)
]
}
,
t
=
2,3
,
…
,
r
f
N
(
t
)
=
b
2
[
X
2
(
1
)
(
t
)
]
β
2
+
b
3
[
X
3
(
1
)
(
t
)
]
β
3
+
⋯
+
b
n
[
X
n
(
1
)
(
t
)
]
β
n
+
u
,
t
=
1,2
,
…
,
r
(
b
1
,
b
2
,
…
,
b
n
,
u
)
T
=
(
B
T
B
)
-
1
B
T
Y
RThe optimisation problem above can be solved by proprietary software. Onceβ
i
=
1, i
=
2,3
,
…
,
n is confirmed, the structure parameters b
1
,
b
2
,
…
,
b
n, u of the model are also determined. Then, the discrete function f
N
(
t
) can be derived by back-substituting the derived results into (26): final predictions can then be obtained by applying f
N
(
t
) to (35) to (37).


## 3. Forecasting China’s Industrial SO2 Emissions

Industry is the dominant factor in promoting China’s economic development. However, SO2 emissions discharged during increasing industrialisation pose a significant threat to the eco-system and biosphere. The effective management of industrial SO2 emissions demands accurate forecasts thereof. The paper argues that industrial output activity levels are considered to be the main factor affecting industrial SO2 emissions. Local governments and enterprises, to increase industrial output, continue to expand the scale of production scale and thereby increase SO2 emissions. China’s industrial SO2 emissions and gross industrial output levels from 2003 to 2010 are given in Table 1: Figure 1 shows the trends in China’s industrial SO2 emissions and gross industrial output. As seen in Figure 1, China’s industrial SO2 emissions and industrial output values are not simply linearly related. With the increasing industrial output, SO2 emissions initially increased and then decreased. Generally, the reasons for this phenomenon are complex; multifarious and interrelated but basically may be summarised as combinations of the following: structural economic change, technical progress, different modalities of demand, and more effective regulatory regimes, and so forth.Table 1
China’s industrial SO2 emissions and gross industrial output values of 2003–2010.
Time (year)
1 (2003)
2 (2004)
3 (2005)
4 (2006)
5 (2007)
6 (2008)
7 (2009)
8 (2010)
Industrial SO2  emissions (104 Tons)
2158.50
2254.90
2549.40
2588.80
2468.09
2321.23
2214.40
2185.15
Industrial output (RMB¥100 Million)
54945.53
65210.03
77230.78
91310.94
110534.88
130260.24
135239.95
160867.01
N.B. Data sourced from the Chinese National Bureau of Statistics Website (http://www.stats.gov.cn/).Figure 1
China’s industrial SO2 emissions correlated with gross industrial output values.GMC (1, 2) and NGMC (1, 2) models were used to forecast—from known industrial output values—China’s SO2 emissions as a time series. Through comparing the modelling and prediction accuracy of the two models, the effectiveness of the proposed NGMC (1, 2) was demonstrated.


### 3.1. Forecasting by GMC (1, 2)

Applying to the GMC (1, 2) model of (7) to (18), the values of parameters n, r, r
f, and r
p in (7), estimates of model parameters in (10) can be obtained and are listed in Table 2. The GMC (1, 2) model from (7) becomes
(41)
d
X
1
(
1
)
(
t
)
d
t
-
0.557995
X
1
(
1
)
(
t
)
=
-
0.015596
X
2
(
1
)
(
t
)
+
1806.203
,
t
=
1,2
,
…
,
8
.Table 2
The values of parametersn
,
r, r
f, and r
p in (7), the estimates of model parameters b
1, b
2, and u in (10), and the values of RMSPEPR and RMSPEPO in (38) and (39).
Parameters
Values
n
2
r
5
r
f
3
r
p
0
b
1
−0.557995
b
2
−0.015596
u
1806.203
RMSPEPR (%)
8.02
RMSPEPO (%)
185.25In summary, the right-hand side of (7), the discrete function f
(
t
) in (13) for the GMC (1, 2) model is obtained and listed in Table 3; the values of RMSPEPR and RMSPEPO in (38) and (39), respectively, are also listed in Table 2. The modelling values and forecasts for China’s industrial SO2 emissions by GMC (1, 2) are listed in Table 4.Table 3
The discrete functionf
(
t
) in (13).
t
f
(
t
)
1
949.250
2
−67.792
3
−1272.31
4
−2696.44
5
−4420.38
6
−6451.98
7
−8561.23
8
−11070.18Table 4
The modelling values and forecasts of China’s industrial SO2 emissions by GMC (1, 2).
Time
Modelling values or forecast
Observations
Percent errors (%)
1 (2003)
2158.50
2158.50
0.00
2 (2004)
2195.32
2254.90
−2.64
3 (2005)
2367.35
2549.40
−7.14
4 (2006)
2398.87
2588.80
−7.34
5 (2007)
2110.67
2468.09
−14.48
6 (2008)
1205.64
2321.23
−48.06
7 (2009)
−630.24
2214.40
−128.46
8 (2010)
−4153.31
2185.15
−290.07As seen in Table2, the root mean squared percentage error for the presample period (RMSPEPR) was 8.02%, but the root mean squared percentage error for the postsample period (RMSPEPO) was as high as 185%. This result was unacceptable. Table 4 shows that the relative errors of modelling and prediction gradually increased over time, especially the relative errors for the period 2008 to 2010, being up to −48%, −128%, and −290%, respectively. This indicated that the traditional GMC (1, 2) model cannot be used to describe the relationship between China’s industrial SO2 emissions and industrial output values nor could it be used to predict future industrial SO2 emissions.


### 3.2. Forecasting by NGMC (1, 2)

Applying the NGMC (1, 2) model of (19) to (37), the values of parameters n
,
r, r
f, and r
p in (21), the estimates of model parameters b
1, b
2 and u in (23), the optimised parameter β
2 in (40) can be obtained and values are listed in Table 5. The NGMC (1, 2) model from (21) becomes
(42)
d
X
1
(
1
)
(
t
)
d
t
-
0.222770
X
1
(
1
)
(
t
)
=
-
46938.3
(
X
2
(
1
)
(
t
)
)
-
0.06305
+
25887.98
,
t
=
1,2
,
…
,
8
.Table 5
The values of parametersn
,
r, r
f, and r
p in (21), the estimates of model parameters b
1, b
2, u in  (23) and the optimized parameter β
2 in (40), and  the values of RMSPEPR and RMSPEPO in (38) and (39).
Parameters
Values
n
2
r
5
r
f
3
r
p
0
b
1
−0.222770
b
2
−46938.3
u
25887.98
β
2
−0.06305
RMSPEPR (%)
2.44
RMSPEPO (%)
5.48In summary, the right-hand side of (21), the discrete function f
N
(
t
) in (26) for the NGMC (
1,2
) model, is obtained and values are listed in Table 6; the values of RMSPEPR and RMSPEPO in (38) and (39), respectively, are also listed in Table 5. The modelling values and forecasts of China’s industrial SO2 emissions by NGMC (
1
,
2
) are listed in Table 7. Among them, the values for 2003 to 2007 are derived by modelling and those for 2008 to 2010 are forecasts.Table 6
The discrete functionf
N
(
t
) in (26).
t
f
N
(
t
)
1
2301.635
2
3437.039
3
4128.824
4
4644.269
5
5074.072
6
5441.376
7
5732.536
8
6006.089Table 7
The modelling values and forecasts of China’s industrial SO2 emissions by NGMC (1,2).
Time
Modelling values or forecast
Observations
Percent errors (%)
1 (2003)
2158.50
2158.50
0.00
2 (2004)
2135.84
2254.90
−5.28
3 (2005)
2526.60
2549.40
−0.89
4 (2006)
2562.03
2588.80
−1.03
5 (2007)
2473.20
2468.09
0.21
6 (2008)
2335.84
2321.23
0.63
7 (2009)
2163.90
2214.40
−2.28
8 (2010)
1984.36
2185.15
−9.19As seen in Table5, the RMSPEPR and RMSPEPO of NGMC (1, 2) were 2.44% and 5.48%, respectively, significantly less than the corresponding RMSPEPR and RMSPEPO values for the GMC (1, 2) model as shown in Table 2. Table 7 also shows the modelling and prediction values for the period 2003 to 2010 by NGMC (1, 2) to be close to the actual values. This indicated that the traditional NGMC (1, 2) model could effectively describe the relationship between China’s industrial SO2 emissions and industrial output values and accurately forecast industrial SO2 emissions.


## 3.1. Forecasting by GMC (1, 2)

Applying to the GMC (1, 2) model of (7) to (18), the values of parameters n, r, r
f, and r
p in (7), estimates of model parameters in (10) can be obtained and are listed in Table 2. The GMC (1, 2) model from (7) becomes
(41)
d
X
1
(
1
)
(
t
)
d
t
-
0.557995
X
1
(
1
)
(
t
)
=
-
0.015596
X
2
(
1
)
(
t
)
+
1806.203
,
t
=
1,2
,
…
,
8
.Table 2
The values of parametersn
,
r, r
f, and r
p in (7), the estimates of model parameters b
1, b
2, and u in (10), and the values of RMSPEPR and RMSPEPO in (38) and (39).
Parameters
Values
n
2
r
5
r
f
3
r
p
0
b
1
−0.557995
b
2
−0.015596
u
1806.203
RMSPEPR (%)
8.02
RMSPEPO (%)
185.25In summary, the right-hand side of (7), the discrete function f
(
t
) in (13) for the GMC (1, 2) model is obtained and listed in Table 3; the values of RMSPEPR and RMSPEPO in (38) and (39), respectively, are also listed in Table 2. The modelling values and forecasts for China’s industrial SO2 emissions by GMC (1, 2) are listed in Table 4.Table 3
The discrete functionf
(
t
) in (13).
t
f
(
t
)
1
949.250
2
−67.792
3
−1272.31
4
−2696.44
5
−4420.38
6
−6451.98
7
−8561.23
8
−11070.18Table 4
The modelling values and forecasts of China’s industrial SO2 emissions by GMC (1, 2).
Time
Modelling values or forecast
Observations
Percent errors (%)
1 (2003)
2158.50
2158.50
0.00
2 (2004)
2195.32
2254.90
−2.64
3 (2005)
2367.35
2549.40
−7.14
4 (2006)
2398.87
2588.80
−7.34
5 (2007)
2110.67
2468.09
−14.48
6 (2008)
1205.64
2321.23
−48.06
7 (2009)
−630.24
2214.40
−128.46
8 (2010)
−4153.31
2185.15
−290.07As seen in Table2, the root mean squared percentage error for the presample period (RMSPEPR) was 8.02%, but the root mean squared percentage error for the postsample period (RMSPEPO) was as high as 185%. This result was unacceptable. Table 4 shows that the relative errors of modelling and prediction gradually increased over time, especially the relative errors for the period 2008 to 2010, being up to −48%, −128%, and −290%, respectively. This indicated that the traditional GMC (1, 2) model cannot be used to describe the relationship between China’s industrial SO2 emissions and industrial output values nor could it be used to predict future industrial SO2 emissions.


## 3.2. Forecasting by NGMC (1, 2)

Applying the NGMC (1, 2) model of (19) to (37), the values of parameters n
,
r, r
f, and r
p in (21), the estimates of model parameters b
1, b
2 and u in (23), the optimised parameter β
2 in (40) can be obtained and values are listed in Table 5. The NGMC (1, 2) model from (21) becomes
(42)
d
X
1
(
1
)
(
t
)
d
t
-
0.222770
X
1
(
1
)
(
t
)
=
-
46938.3
(
X
2
(
1
)
(
t
)
)
-
0.06305
+
25887.98
,
t
=
1,2
,
…
,
8
.Table 5
The values of parametersn
,
r, r
f, and r
p in (21), the estimates of model parameters b
1, b
2, u in  (23) and the optimized parameter β
2 in (40), and  the values of RMSPEPR and RMSPEPO in (38) and (39).
Parameters
Values
n
2
r
5
r
f
3
r
p
0
b
1
−0.222770
b
2
−46938.3
u
25887.98
β
2
−0.06305
RMSPEPR (%)
2.44
RMSPEPO (%)
5.48In summary, the right-hand side of (21), the discrete function f
N
(
t
) in (26) for the NGMC (
1,2
) model, is obtained and values are listed in Table 6; the values of RMSPEPR and RMSPEPO in (38) and (39), respectively, are also listed in Table 5. The modelling values and forecasts of China’s industrial SO2 emissions by NGMC (
1
,
2
) are listed in Table 7. Among them, the values for 2003 to 2007 are derived by modelling and those for 2008 to 2010 are forecasts.Table 6
The discrete functionf
N
(
t
) in (26).
t
f
N
(
t
)
1
2301.635
2
3437.039
3
4128.824
4
4644.269
5
5074.072
6
5441.376
7
5732.536
8
6006.089Table 7
The modelling values and forecasts of China’s industrial SO2 emissions by NGMC (1,2).
Time
Modelling values or forecast
Observations
Percent errors (%)
1 (2003)
2158.50
2158.50
0.00
2 (2004)
2135.84
2254.90
−5.28
3 (2005)
2526.60
2549.40
−0.89
4 (2006)
2562.03
2588.80
−1.03
5 (2007)
2473.20
2468.09
0.21
6 (2008)
2335.84
2321.23
0.63
7 (2009)
2163.90
2214.40
−2.28
8 (2010)
1984.36
2185.15
−9.19As seen in Table5, the RMSPEPR and RMSPEPO of NGMC (1, 2) were 2.44% and 5.48%, respectively, significantly less than the corresponding RMSPEPR and RMSPEPO values for the GMC (1, 2) model as shown in Table 2. Table 7 also shows the modelling and prediction values for the period 2003 to 2010 by NGMC (1, 2) to be close to the actual values. This indicated that the traditional NGMC (1, 2) model could effectively describe the relationship between China’s industrial SO2 emissions and industrial output values and accurately forecast industrial SO2 emissions.


## 4. Conclusions

Grey prediction modelling with the addition of a convolution integral to GMC(
1
,
n
) system greatly improved the prediction accuracy of the traditional multiply-grey model GM (
1
,
n
) and was successfully applied in practice. The current GMC (
1
,
n
) model is linear; therefore, the premise of obtaining high-precision modelling results by GMC (
1
,
n
) is the existence of similar trends between the predicted variables and their influencing factors. The research introduced a power exponent into a 1-AGO series of relative factors of GMC (
1
,
n
) to reflect the nonlinear interaction of the associated series with that predicted. The nonlinear optimisation model is built to solve location parameters. The prediction examples related to China’s industrial SO2 emissions showed that the GMC (1, 2) model cannot effectively describe the nonlinear relationship between China’s industrial SO2 emissions and their industrial output: the ensuing prediction errors were insurmountable. However, the NGMC (1, 2) model proposed in this paper effectively described the relationship and achieved satisfactory prediction accuracy.Compared with the traditional GMC(
1
,
n
),  n
-
1power exponents of the predicted variables’ relevant factors β
2
,
β
3
,
…
,
β
n are introduced into the NGMC (
1
,
n
) model to reflect the effect of these upon nonlinear system behaviours and interactions. The unknown parameters are determined by a computer program, which calculates the minimum average relative percentage error of the forecasting model. This strengthens the adaptability of the NGMC (
1
,
n
) model towards the original data and eventually improves the forecast accuracy.



---

*Source: 580161-2014-03-24.xml*