*Article* **Stability of Time Series Models Based on Fractional-Order Weakening Buffer Operators**

**Chong Li 1, Yingjie Yang <sup>2</sup> and Xinping Zhu 1,\***


**Abstract:** Different weakening buffer operators in a time-series model analysis usually result in different model sensitivities, which sometimes affect the effectiveness of relevant operator-based methods. In this paper, the stability of two classic fractional-order weakening buffer operator-based series models is studied; then, a new data preprocessing method based on a novel fractional-order bidirectional weakening buffer operator is provided, whose effect in improving the model's stability is tested and utilized in prediction problems. Practical examples are employed to demonstrate the efficiency of the proposed method in improving the model's stability in noise scenarios. The comparison indicates that the proposed method overcomes the disadvantage of many weakening buffer operators in the subjectively biased weighting of the new or old information in forecasting. These expand the application of the proposed method in time series analysis.

**Keywords:** time series; sequence operator; fractional-order operator; model stability; model perturbation analysis

### **1. Introduction**

Time series prediction models play important roles in many real-world applications, including economic and social forecasting problems [1,2]. A literature survey reveals that fractional-order dynamical models sometimes provide more valuable insights to portray complex natural phenomena than their integer-order counterparts [3,4]. However, data noise and missing data are often encountered in model applications, and they affect the effectiveness of prediction results. The existence and stability of stationary solutions are basic problems for models subject to exogenous random perturbations [5,6]. Researchers have proposed ways and methods to improve the accuracy of model solutions [7–10].

Sequence operators play an important role in reducing the interference information in collected data and highlight the potential development process of analyzed objects. Among them are the widely used exponential smoothing operators [11], especially the seasonal exponential smoothing method in the application of forecasting [12,13]. The buffer operators in the grey systems theory also exhibit similar features [14]. They are based on the three axioms of the buffer operator and are suitable for small sample analyses. They greatly improve the performance of grey prediction models in real applications [15]. These operators can be further subdivided into weakening operators and strengthening operators [16,17]. When it comes to system development trend analysis or series prediction, weakening buffer operators (WBO) are preferred. These types of operators include the average weakening buffer operator (AWBO), the geometric averageweakening buffer operator (GAWBO), and the weighted WBOs (WAWBO, WGAWBO), etc., see [18,19]. The main disadvantage of these operators is their subjective determination of the weight coefficients of data [20]. This will undoubtedly miss some of the useful information under certain conditions and limit their applications.

Wu et al. first studied the essence of WBOs in grey prediction models [21–23]. Based on the perturbation theory, they found that series prediction models are more sensitive

**Citation:** Li, C.; Yang, Y.; Zhu, X. Stability of Time Series Models Based on Fractional-Order Weakening Buffer Operators. *Fractal Fract.* **2023**, *7*, 554. https://doi.org/10.3390/ fractalfract7070554

Academic Editors: Angelo B. Mingarelli, Leila Gholizadeh Zivlaei and Mohammad Dehghan

Received: 17 May 2023 Revised: 26 June 2023 Accepted: 27 June 2023 Published: 17 July 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to earlier data than newer data in a given sequence. Their findings reveal that the classic integral accumulated generating operator, 1-AGO, and the extended fractional accumulated generating operator attach more importance to the earlier data than to recent data, and both operators donot satisfy the priority theory of new information in the grey systems theory. To solve this problem, they later proposed some reverse fractional-order operators in Wu et al. [23] and Wu et al. [21], which attach more importance to the morerecent data and thus are consistent with the priority theory of new information. Recent studies show that these new fractional-order accumulation and inverse accumulation operators could improve the prediction performance of time series models [24–26]. However, these new fractional-order operators emphasize either the new or old data and thus are only suitable for series analyses with special time preferences.

To solve the above problems, this paper applies a data preprocessing method based on a novel fractional-order bidirectional weakening buffer operator, which was first proposed by Li et al. [27]. Recent studies have demonstrated the good performance of this fractionalorder weakening buffer operator in predictions against data noise interference [28]. The effect of this adopted operator, together with two other widely used fractional-order operators, on improving the stability of time series models is examined and compared. When applying the perturbation theory, this study shows that: (1) the disadvantage of the two classic fractional-order weakening buffer operators lies in the subjectively biased weighting of either the new or old information in forecasting; (2) the proposed data preprocessing method, based on a fractional bidirectional weakening buffer operator, improves the stability of time series models; (3) this newly adopted fractional-order weakening buffer operator-based method can more objectively deal with both the old and new data noise in samples and thereby improve the accuracy of model predictions. Both theoretical investigations and numerical experiments show favorable effects for the new fractionalorder bidirectional weakening buffer operator in improving the performance of time series models. These enlarge the application of series prediction models.

This paper is structured as follows: Section 2 introduces three important fractionalorder weakening buffer operators used in the series models. Section 3 analyzes the performance of the proposed operator-based data preprocessing method on the sensitivity of a fractional-order accumulate discrete grey model *DGMp*(1, 1). Section 4 shows the model perturbation analysis result of the proposed method and its comparison operator on another model—the new information prior to the grey model *NIGM*(*p*, 1). The case studiesin Section 5 demonstrate the effect of the new method in improving the model's stability. Finally, conclusions are drawn in Section 6.

#### **2. Fractional-Order Weakening Buffer Operators**

Weakening buffer operators are commonly used in series prediction models for finding the implicit pattern in samples. In order to show the advantage of the adopted fractionalorder bidirectional weakening buffer operator in the proposed data preprocessing method, two other widely used classic fractional-order weakening buffer operators are introduced in Definition 1 and Definition 2. They are chosen as the comparative objects of the proposed one.

**Definition 1** [22]**.** *Let X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / *be the original sequence, and <sup>D</sup>pX*(0) <sup>=</sup> *<sup>X</sup><sup>p</sup>* <sup>=</sup> {*xp*(1), *<sup>x</sup>p*(2),..., *<sup>x</sup>p*(*n*)} *be the fractional <sup>p</sup>*-order (<sup>0</sup> <sup>&</sup>lt; *<sup>p</sup>* <sup>&</sup>lt; <sup>1</sup>) *accumulated generating operator sequence of non-negative X*(0)*, with operator D<sup>p</sup> satisfying:*

$$\mathbf{x}^{p}(k) = \sum\_{i=1}^{k} \mathbb{C}\_{k-i+p-1}^{k-i} \mathbf{x}^{(0)}(i), \qquad (k = 1, 2, \dots, n) \tag{1}$$

*where C*<sup>0</sup> *<sup>p</sup>*−<sup>1</sup> <sup>=</sup> 1, *<sup>C</sup><sup>k</sup> <sup>k</sup>*−<sup>1</sup> <sup>=</sup> 0, *<sup>C</sup>k*−*<sup>i</sup> <sup>k</sup>*−*i*+*p*−<sup>1</sup> <sup>=</sup> (*k*−*i*+*p*−1)(*k*−*i*+*p*−2)...(*p*+1)(*p*) (*k*−*i*)! . **Definition 2** [23]**.** *Let X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / *be the original sequence, and DrX*(0) = *X*(*r*) = . *x*(*r*)(1), *x*(*r*)(2),..., *x*(*r*)(*n*) / *be defined as the reverse fractional r*-order *accumulated generating operator sequence of non-negative X*(0)*, with operator D<sup>r</sup> satisfying:*

$$\mathbf{x}\_{(r)}(k) = \sum\_{i=k}^{n} \mathbb{C}\_{i-k+r-1}^{i-k} \mathbf{x}^{(0)}(i), \qquad (k = 1, 2, \dots n) \tag{2}$$

*where C*<sup>0</sup> *<sup>r</sup>*−1<sup>=</sup> 1, *<sup>C</sup><sup>n</sup> <sup>n</sup>*−1<sup>=</sup> 1, *<sup>C</sup>i*−*<sup>k</sup> <sup>i</sup>*−*k*+*r*−<sup>1</sup> <sup>=</sup> (*i*−*k*+*r*−1)(*i*−*k*+*r*−2)···(*r*+1)*<sup>r</sup>* (*i*−*k*)! .

If order parameters *p* in Definition 1 and *r* in Definition 2 are limited to positive integers, the p-AGO and the r-IAGO are obtained. The important properties of these operators have been intensively studied [23,29]. In this study, the fractional-order bidirectional weakening buffer operator adopted in the proposed data preprocessing method of the time series prediction models is provided in the following definition:

**Definition 3** [27]**.** *For the original sequence X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / *, an operator sequence Xv* <sup>=</sup> {*xv*(1), *xv*(2),..., *xv*(*n*)} *is obtained by DvX*(0) <sup>=</sup> *Xv, where the time series operator Dv*(*<sup>v</sup>* <sup>∈</sup> *<sup>R</sup>*+) *is a fractional-order bidirectional weakening buffer operator, and elements in sequence Xv satisfy:*

$$\mathbf{x}\_{\upsilon}(i) = \left(\sum\_{j=i-a(i)}^{i} w(i-j)\mathbf{x}^{(0)}(j) + \sum\_{j=i+1}^{i+a(i)} w(j-i)\mathbf{x}^{(0)}(j)\right) / \left(\sum\_{j=i-a(i)}^{i} w(i-j) + \sum\_{j=i+1}^{i+a(i)} w(j-i)\right),\tag{3}$$

*where α*(*i*) = min(*i* − 1, *n* − *i*)*, and*

$$w(i) = \begin{cases} \left(n^v(v+1-n) + (n-1)^{v+1}\right) / \Gamma(v+2), & i = n\\ \left((i+1)^{v+1} - 2i^{v+1} + (i-1)^{v+1}\right) / \Gamma(v+2), & i = 1, 2, \dots, (n-1) \\\ 1/\Gamma(v+2), & i = 0 \end{cases} \quad i = 1, 2, \dots, (n-1) \quad \text{(4)}$$

It has been proved by Li et al. [27] that the fractional-order bidirectional weakening buffer operator *Dv* defined in (3) is a weakening buffer operator. According to Definition 3, the value adjustment of an element in a given series is determined by both the older and newer data, while it is exclusively determined by the old or new data if applying Definition 1 or Definition 2. In the following sections, we will test the performance of operator *Dv* in improving the stability of two grey prediction models, which are based on the operators presented in Definition 1 and Definition 2, respectively.

#### **3. Perturbation Analysis of Model 1: The Fractional-Order Accumulate Discrete Grey Model in [22]**

Whether a time series model is sensitive to the data samples plays different functions in different applications. In this section, the perturbation theory is employed to test the effect of the fractional-order bidirectional weakening buffer operator *Dv* on the sensitivity of the fractional-order accumulate discrete grey model *DGMp*(1, 1) designed in [22]. In this paper, we denote the consistent matrix norm by . and denote the spectral matrix norm by .2.

**Lemma 1** [30]**.** *Let matrices <sup>A</sup>* <sup>∈</sup> *<sup>C</sup>n*×*n*, *<sup>F</sup>* <sup>∈</sup> *<sup>C</sup>n*×*n, and vectors <sup>b</sup>* <sup>∈</sup> *<sup>C</sup>n*, *<sup>c</sup>* <sup>∈</sup> *<sup>C</sup>n. Assume that thematrix or vector norm* . *is consistent, rank*(*A*) = *rank*(*A* + *F*) = *n and the matrix norm*  *A*−<sup>1</sup> *F* < 1 *is true, then the least squares estimation solution x of the linear system models Ax* = *b and the least squares estimation solution x* + *h of the model* (*A* + *F*)(*x* + *h*) = *b* + *c satisfy:*

$$||h|| \le \frac{\kappa\_{\mathsf{H}}}{\gamma\_{\mathsf{H}}} \left( \frac{||F||\_{2}||\mathsf{x}||}{||A||} + \frac{||c||}{||A||} + \frac{\kappa\_{\mathsf{H}}}{\gamma\_{\mathsf{H}}} \frac{||F||\_{2}||\gamma\_{\mathsf{x}}||}{||A||^{2}} \right) \tag{5}$$

*where κ*† = *<sup>A</sup>*−<sup>1</sup> 2 *A* , *<sup>γ</sup>*† <sup>=</sup> <sup>1</sup> − *A*−1<sup>2</sup> *F* <sup>2</sup>, *γ<sup>x</sup>* = *b* − *Ax*.

To keep the consistency of the parameters, the fractional-order parameter *p*/*q* used in [22] is replaced by parameter *p* in the rest of this paper.

**Definition 4** [22]**.** *Let X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / *be the original sequence, and X<sup>p</sup>* = {*xp*(1), *<sup>x</sup>p*(2),..., *<sup>x</sup>p*(*n*)} *be its fractional p-order accumulated generating operator sequence based on Definition 1, then the following time series model is called a fractional-order accumulate discrete grey model DGMp*(1, 1)*:*

$$\mathbf{x}^{p}(k+1) = \beta\_1 \mathbf{x}^{p}(k) + \beta\_2, \qquad k = 1, 2, \dots, n-1 \tag{6}$$

The least squares estimation of model parameters *β*<sup>1</sup> and *β*<sup>2</sup> in model (6) are:

$$
\begin{bmatrix} \beta\_2\\ \beta\_1 \end{bmatrix} = \left(B^T B\right)^{-1} B^T Y\_\prime \tag{7}
$$

.

where

$$B = \begin{bmatrix} 1 & \mathfrak{x}^p(1) \\ 1 & \mathfrak{x}^p(2) \\ \vdots & \vdots \\ 1 & \mathfrak{x}^p(n-2) \\ 1 & \mathfrak{x}^p(n-1) \end{bmatrix}, \qquad Y = \begin{bmatrix} \mathfrak{x}^p(2) \\ \mathfrak{x}^p(3) \\ \vdots \\ \mathfrak{x}^p(n-1) \\ \mathfrak{x}^p(n) \end{bmatrix}$$

For simplicity, let *β* = [*β*2, *β*1] *<sup>T</sup>*. The perturbation analysis result of model (6) in [22] is summarized in the following Theorems 1–3.

**Theorem 1** [22]**.** *For the DGMp*(1, 1) *model with original data X*(0) = . *x*(0)(1),..., *x*(0)(*n*) / *, let β be its least squares estimation. When the raw data is disturbed, DGMp*(1, 1) *becomes* (*B* + Δ*B*)*β* 1 = (*Y* + Δ*Y*)*, where* Δ*B and* Δ*Y are determined by the disturbance item ε. Assume that the matrixnorm condition <sup>B</sup>*−<sup>1</sup> 2 Δ*B* <sup>2</sup> < 1 *holds for DGMp*(1, 1)*. Let κ*† = *<sup>B</sup>*−<sup>1</sup> 2 *B , γ*† = 1− *<sup>B</sup>*−<sup>1</sup> 2 Δ*B* <sup>2</sup>*, γβ* = *Y* − *Bβ. If the disturbance occurs in the first data of the original sequence, that is <sup>x</sup>*1(0)(1) = *<sup>x</sup>*(0)(1) + *<sup>ε</sup>, then the difference between solutions <sup>β</sup>* 1 *and β is denoted by δ*(*x*(0)(1)) *and it satisfies:*

$$\|\delta(\mathbf{x}^{(0)}(1))\| \le |\varepsilon| \frac{\kappa\_{\mathbf{r}}}{\gamma\_{\mathbf{r}}} \left( \frac{\sqrt{\sum\_{k=1}^{n-1} \left(\mathbf{C}\_{k+p-2}^{k-1}\right)^2} \|\beta\|}{\|B\|} + \frac{\sqrt{\sum\_{k=2}^{n} \left(\mathbf{C}\_{k+p-2}^{k-1}\right)^2}}{\|B\|} + \frac{\kappa\_{\mathbf{r}}}{\gamma\_{\mathbf{r}}} \frac{\sqrt{\sum\_{k=1}^{n-1} \left(\mathbf{C}\_{k+p-2}^{k-1}\right)^2} \|\gamma\_{\beta}\|}{\|B\|} \right). \tag{8}$$

**Theorem 2** [22]**.** *For the same model and parameters as those in Theorem 1, if the disturbance occurs in the non-boundary nodes, that is <sup>x</sup>*1(0)(*k*) = *<sup>x</sup>*(0)(*k*) + *<sup>ε</sup>,* (*<sup>k</sup>* <sup>=</sup> 2, 3, ... , *<sup>n</sup>* <sup>−</sup> <sup>1</sup>)*, then the difference between solutions of model DGMp*(1, 1) *is denoted by δ*(*x*(0)(*k*))*, and it satisfies:*

$$\left\lVert\delta(\mathbf{x}^{(0)}(k))\right\rVert \leq |\mathfrak{e}| \frac{\mathsf{x}\_{\mathsf{t}}}{\gamma\_{\mathsf{t}}} \left( \frac{\sqrt{\sum\_{i=1}^{n-k} \left(\mathsf{C}\_{i+p-2}^{i-1}\right)^{2}} \left\lVert\beta\right\rVert}{\left\lVert\boldsymbol{B}\right\rVert} + \frac{\sqrt{\sum\_{i=1}^{n-k+1} \left(\mathsf{C}\_{i+p-2}^{i-1}\right)^{2}}}{\left\lVert\boldsymbol{B}\right\rVert} + \frac{\mathsf{x}\_{\mathsf{t}}}{\gamma\_{\mathsf{t}}} \frac{\sqrt{\sum\_{i=1}^{n-k} \left(\mathsf{C}\_{i+p-2}^{i-1}\right)^{2}} \left\lVert\gamma\_{\beta}\right\rVert}{\left\lVert\boldsymbol{B}\right\rVert^{2}} \right). \tag{9}$$

**Theorem 3** [22]**.** *For the same model and parameters as those in Theorem 1, if the disturbance occurs in the last data of the original sequence, that is <sup>x</sup>*1(0)(*n*) = *<sup>x</sup>*(0)(*n*) + *<sup>ε</sup>, then the difference between solutions of model DGMp*(1, 1) *is denoted by δ*(*x*(0)(*n*))*, and it satisfies:*

$$\|\delta(\mathbf{x}^{(0)}(n))\| \le \frac{\kappa\_{\!\!\!+}}{\gamma\_{\!\!\!\!+}} \frac{|\varepsilon|}{||B||} \tag{10}$$

Based on the above theorems, Wu et al. [22] demonstrate that the proposed *DGMp*(1, 1) has better solution stability than the classic *GM*(1, 1) model. Now, to test the advantageof the fractional-order bidirectional weakening buffer operator *Dv* in improving the model's stability, we first apply a preprocessing step to the original time series based on operator *Dv* and then carry out the *DGMp*(1, 1) analysis. The perturbation theory is employed to test the effect of operator *Dv* on the solution stability of *DGMp*(1, 1).

Let *X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / be the original sequence, *Xv* = {*xv*(1), *xv*(2), ... , *xv*(*n*)} be its weakening buffer operator sequence based on operator *Dv*. Apply the *DGMp*(1, 1) model to sequence *Xv* and let *β<sup>v</sup>* be the new model solution values of parameters [*β*2, *β*1] *<sup>T</sup>*; then we have

$$\beta\_{\upsilon} = \left(B\_{\upsilon}^{\,T} B\_{\upsilon}\right)^{-1} B\_{\upsilon}^{\,T} Y\_{\upsilon} \,. \tag{11}$$

where *x<sup>p</sup> <sup>v</sup>* (*k*) = ∑*<sup>k</sup> <sup>i</sup>*=<sup>1</sup> *<sup>C</sup>k*−*<sup>i</sup> <sup>k</sup>*−*i*+*p*−1*xv*(*i*), and

$$B\_{\upsilon} = \begin{bmatrix} 1 & \boldsymbol{\varkappa}\_{\upsilon}^{p}(1) \\ 1 & \boldsymbol{\varkappa}\_{\upsilon}^{p}(2) \\ \vdots & \vdots \\ 1 & \boldsymbol{\varkappa}\_{\upsilon}^{p}(n-2) \\ 1 & \boldsymbol{\varkappa}\_{\upsilon}^{p}(n-1) \end{bmatrix}, \qquad Y\_{\upsilon} = \begin{bmatrix} \boldsymbol{\varkappa}\_{\upsilon}^{p}(2) \\ \boldsymbol{\varkappa}\_{\upsilon}^{p}(3) \\ \vdots \\ \boldsymbol{\varkappa}\_{\upsilon}^{p}(n-1) \\ \boldsymbol{\varkappa}\_{\upsilon}^{p}(n) \end{bmatrix}.$$

When the raw data is disturbed, let Δ*Bv* and Δ*Yv* be the disturbance items of *Bv* and *Yv*, respectively, and *β* <sup>1</sup>*<sup>v</sup>* be the new model solution. Assume that *<sup>B</sup>*−<sup>1</sup> *<sup>v</sup>* 2 Δ*Bv* <sup>2</sup> < 1 and set *κv*† = *<sup>B</sup>*−<sup>1</sup> *<sup>v</sup>* 2 *Bv* , *γv*† = 1− *<sup>B</sup>*−<sup>1</sup> *<sup>v</sup>* 2 Δ*Bv* <sup>2</sup> and *γv<sup>β</sup>* = *Yv* − *Bvβv*. The parameters *α*(*i*) and *w*(*i*) used in the remaining part of this section are the same as those defined in model (3). The main results are stated in the following theorems.

**Theorem 4.** *If a raw data disturbance occurs in the first data of the original sequence, that is <sup>x</sup>*1(0)(1) = *<sup>x</sup>*(0)(1) + *<sup>ε</sup>, then the difference between <sup>β</sup>* 1*<sup>v</sup> and β<sup>v</sup> satisfies:*

$$\left\|\delta(\mathbf{x}^{(0)}(1))\right\| \leq \left|\varepsilon\right| \frac{\kappa\_{\upsilon\mathsf{t}}}{\gamma\_{\upsilon\mathsf{t}}} \left(\frac{M\_{\upsilon}||\beta\_{\upsilon}||}{||B\_{\upsilon}||} + \frac{N\_{\upsilon}}{||B\_{\upsilon}||} + \frac{\kappa\_{\upsilon\mathsf{t}}}{\gamma\_{\upsilon\mathsf{t}}} \frac{M\_{\upsilon}||\gamma\_{\upsilon\beta}||}{||B\_{\upsilon}||^{2}}\right),\tag{12}$$

*where: (i). If sequence size n is even, then*:

$$\begin{array}{l} M\_{\upsilon} = \sqrt{\frac{\sum\_{i=1}^{\upsilon/2} \left( \sum\_{j=1}^{i} \mathcal{C}\_{j+p-2}^{j-1} \frac{w(i-j)}{w(0) + \sum\_{k=1}^{a(i-j+1)} w(k)} \right)^{2} + \sum\_{i=n/2+1}^{n-1} \left( \sum\_{j=1}^{n/2} \mathcal{C}\_{i+p-j-1}^{i-j} \frac{w(j-1)}{w(0) + \sum\_{k=1}^{a(j)} w(k)} \right)^{2}}{\sqrt{\sum\_{i=2}^{n/2} \left( \sum\_{j=1}^{i} \mathcal{C}\_{j+p-2}^{j-1} \frac{w(i-j)}{w(0) + 2\sum\_{k=1}^{a(i-j+1)} w(k)} \right)^{2} + \sum\_{i=n/2+1}^{n} \left( \sum\_{j=1}^{n/2} \mathcal{C}\_{i+p-j-1}^{i-j} \frac{w(j-1)}{w(0) + 2\sum\_{k=1}^{a(j)} w(k)} \right)^{2}} \tag{13}$$

*(ii). If sequence size n is odd, then:*

$$\begin{array}{l} \mathcal{M}\_{\mathcal{V}} = \sqrt{\frac{(n+1)/2}{\sum\limits\_{i=1}^{\left(1+\sum\limits\_{j=1}^{i}\frac{\text{w}}{\text{w}}\right)-\sum\limits\_{k=1}^{\left(1-\frac{\text{w}}{\text{w}}\right)-\text{w}}\text{w}(k)}} + \sum\limits\_{i=\left(n+3\right)/2}^{n-1} \left(\sum\limits\_{j=1}^{\left(n+1\right)/2} \mathcal{C}^{i-j}\_{i+p-j-1} \frac{\text{w}(j-1)}{\text{w}(0)+2\sum\limits\_{k=1}^{\left(1\right)}\text{w}(k)}\right)^{2} \\ \mathcal{N}\_{\mathcal{V}} = \sqrt{\frac{(n+1)/2}{\sum\limits\_{i=1}^{\left(n+1\right)}\left(\sum\limits\_{j=1}^{i}\frac{\text{w}(i-j)}{\text{w}(0)+2\sum\limits\_{k=1}^{\left(n-j+1\right)}\text{w}(k)}\right)^{2} + \sum\limits\_{i=\left(n+3\right)/2}^{n} \left(\sum\limits\_{j=1}^{(n+1)/2} \mathcal{C}^{i-j}\_{i+p-j-1} \frac{\text{w}(j-1)}{\text{w}(0)+2\sum\limits\_{k=1}^{\left(n-j\right)}\text{w}(k)}\right)^{2} . \end{array} \tag{14}$$

**Proof of Theorem 4.** (i). Sequence size *<sup>n</sup>* is even. When *<sup>x</sup>*1(0)(1) = *<sup>x</sup>*(0)(1) + *<sup>ε</sup>*, according to the logic of *DGMp*(1, 1), we have

Δ*Bv* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 *ε* 0 *C*<sup>1</sup> *<sup>p</sup><sup>ε</sup>* <sup>+</sup> *<sup>w</sup>*(1)*<sup>ε</sup> w*(0)+2*w*(1) . . . . . . 0 *Cn*/2−<sup>1</sup> *<sup>p</sup>*+*n*/2−2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*/2−<sup>2</sup> *p*+*n*/2−3 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>w</sup>*(*n*/2−1)*<sup>ε</sup> w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) 0 *Cn*/2 *<sup>p</sup>*+*n*/2−1*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*/2−<sup>1</sup> *p*+*n*/2−2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>1</sup> *<sup>p</sup> <sup>w</sup>*(*n*/2−1)*<sup>ε</sup> w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) 0 *Cn*/2+<sup>1</sup> *<sup>p</sup>*+*n*/2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*/2 *p*+*n*/2−1 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>2</sup> *p*+1 *w*(*n*/2−1)*ε w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) . . . . . . <sup>0</sup> *<sup>C</sup>n*−<sup>2</sup> *<sup>p</sup>*+*n*−3*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*−<sup>3</sup> *<sup>p</sup>*+*n*−<sup>4</sup> *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>n*/2−<sup>1</sup> *p*+*n*/2−2 *w*(*n*/2−1)*ε w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Δ*Yv* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *C*1 *<sup>p</sup><sup>ε</sup>* <sup>+</sup> *<sup>w</sup>*(1)*<sup>ε</sup> w*(0)+2*w*(1) *C*2 *<sup>p</sup>*+1*<sup>ε</sup>* + *<sup>C</sup>*<sup>1</sup> *<sup>p</sup> <sup>w</sup>*(1)*<sup>ε</sup> <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> *<sup>w</sup>*(2)*<sup>ε</sup> w*(0)+2(*w*(1)+*w*(2)) . . . *Cn*/2 *<sup>p</sup>*+*n*/2−1*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*/2−<sup>1</sup> *p*+*n*/2−2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>1</sup> *<sup>p</sup> <sup>w</sup>*(*n*/2−1)*<sup>ε</sup> w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) *Cn*/2+<sup>1</sup> *<sup>p</sup>*+*n*/2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*/2 *p*+*n*/2−1 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>2</sup> *p*+1 *w*(*n*/2−1)*ε w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) *Cn*/2+<sup>2</sup> *<sup>p</sup>*+*n*/2+1*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*/2+<sup>1</sup> *p*+*n*/2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>3</sup> *p*+2 *w*(*n*/2−1)*ε w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) . . . *<sup>C</sup>n*−<sup>1</sup> *<sup>p</sup>*+*n*−2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*−<sup>2</sup> *<sup>p</sup>*+*n*−<sup>3</sup> *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>n*/2 *p*+*n*/2−1 *w*(*n*/2−1)*ε w*(0)+2 *α*(*n*/2) ∑ *k*=1 *w*(*k*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

,

Now, calculate the norms of the matrices Δ*Bv* and Δ*Yv*, extract parameter *ε* and apply Lemma 1; parameter expressions in Equation (13) are then obtained.

(ii). When sequence size *n* is odd, we have:

Δ*Bv* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 *ε* 0 *C*<sup>1</sup> *<sup>p</sup><sup>ε</sup>* <sup>+</sup> *<sup>w</sup>*(1)*<sup>ε</sup> w*(0)+2*w*(1) . . . . . . 0 *C*(*n*−1)/2 *<sup>p</sup>*+(*n*−3)/2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>*(*n*−3)/2 *p*+(*n*−5)/2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>w</sup>*((*n*−1)/2)*<sup>ε</sup> w*(0)+2 *α*((*n*+1)/2) ∑ *k*=1 *w*(*k*) 0 *C*(*n*+1)/2 *<sup>p</sup>*+(*n*−1)/2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>*(*n*−1)/2 *p*+(*n*−3)/2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>1</sup> *<sup>p</sup> <sup>w</sup>*((*n*−1)/2)*<sup>ε</sup> w*(0)+2 *α*((*n*+1)/2) ∑ *k*=1 *w*(*k*) . . . . . . <sup>0</sup> *<sup>C</sup>n*−<sup>2</sup> *<sup>p</sup>*+*n*−3*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*−<sup>3</sup> *<sup>p</sup>*+*n*−<sup>4</sup> *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*(*n*−3)/2 *p*+(*n*−5)/2 *w*((*n*−1)/2)*ε w*(0)+2 *α*((*n*+1)/2) ∑ *k*=1 *w*(*k*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Δ*Yv* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *C*1 *<sup>p</sup><sup>ε</sup>* <sup>+</sup> *<sup>w</sup>*(1)*<sup>ε</sup> w*(0)+2*w*(1) *C*2 *<sup>p</sup>*+1*<sup>ε</sup>* + *<sup>C</sup>*<sup>1</sup> *<sup>p</sup> <sup>w</sup>*(1)*<sup>ε</sup> <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> *<sup>w</sup>*(2)*<sup>ε</sup> w*(0)+2*w*(1)+2*w*(2) . . . *C*(*n*+1)/2 *<sup>p</sup>*+(*n*−1)/2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>*(*n*−1)/2 *p*+(*n*−3)/2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>1</sup> *<sup>p</sup> <sup>w</sup>*((*n*−1)/2)*<sup>ε</sup> w*(0)+2 *α*((*n*+1)/2) ∑ *k*=1 *w*(*k*) *C*(*n*+3)/2 *<sup>p</sup>*+(*n*+1)/2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>*(*n*+1)/2 *p*+(*n*−1)/2 *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*<sup>2</sup> *p*+1 *w*((*n*−1)/2)*ε w*(0)+2 *α*((*n*+1)/2) ∑ *k*=1 *w*(*k*) . . . *<sup>C</sup>n*−<sup>1</sup> *<sup>p</sup>*+*n*−2*<sup>ε</sup>* <sup>+</sup> *<sup>C</sup>n*−<sup>2</sup> *<sup>p</sup>*+*n*−<sup>3</sup> *w*(1)*ε <sup>w</sup>*(0)+2*w*(1) <sup>+</sup> ··· *<sup>C</sup>*(*n*−1)/2 *p*+(*n*−3)/2 *w*((*n*−1)/2)*ε w*(0)+2 *α*((*n*+1)/2) ∑ *k*=1 *w*(*k*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

,

Then, calculate the matrix norm ||Δ*Bv*||<sup>2</sup> and vector norm ||Δ*Yv*||, extract parameter *ε* and apply Lemma 1; parameter expressions in Equation (14) are then obtained.

**Theorem 5.** *If a raw data disturbance occurs in the non-boundary nodes of the first half of the original sequence, that is <sup>x</sup>*1(0)(*k*) = *<sup>x</sup>*(0)(*k*) + *<sup>ε</sup>*, *then the difference between <sup>β</sup>* 1*<sup>v</sup> and β<sup>v</sup> satisfies expression (12); however, the values of parameters Mv and Nv are determined by the following scenarios. Let ϕ*<sup>1</sup> *be the largest integer less than k*/2*, ϕ*<sup>2</sup> *be the largest integer less than* (*k* − 1)/2*, then Mv and Nv are:*

*(i). If sequence size n is even, then k* = 2, 3, . . . , *n*/2 *and*

$$\begin{array}{l} M\_{\mathrm{F}} = \sqrt{\frac{\sum\_{i=1}^{2+q} \binom{i-q}{i}}{\sum\_{j=q+1}^{n} \left( \sum\_{j=1}^{i-1} \mathcal{C}\_{j+p-2}^{i-1} \frac{w((k+j-i-1)}{w(0)+2\zeta\_{\mathrm{u}=1}^{i(j-1)}w(k)} \right)^{2} + \sum\_{i=n/2+q+1}^{n-1} \left( \sum\_{j=1}^{n/2} \mathcal{C}\_{p+i+j-q-1-n/2}^{i+j-n/2} \frac{w((k+j-q-1-n/2))}{w(0)+2\zeta\_{\mathrm{u}=1}^{i(n/2+q+1-n/2)} w(k)} \right)^{2}}{\sum\_{i=q+1}^{n/2+q} \left( \sum\_{j=1}^{i-q} \mathcal{C}\_{j+p-2}^{i-1} \frac{w((k+j-i-1)-1)}{w(0)+2\zeta\_{\mathrm{u}=1}^{i(j-n/2)} w(k)} \right)^{2} + \sum\_{i=n/2+q+1}^{n} \left( \sum\_{j=1}^{n/2} \mathcal{C}\_{p+i+j-q-1-n/2}^{i+j-n/2} \frac{w((k+j-q-1-n/2))}{w(0)+2\zeta\_{\mathrm{u}=1}^{i(n/2+q+1-n/2)} w(k)} \right)^{2}} \end{array} \tag{15}$$

*(ii). If sequence size n is odd, then k* = 2, 3, . . . ,(*n* + 1)/2 *and*

$$\begin{array}{l} \mathcal{M}\_{\upsilon} = \sqrt{\frac{(n+1)/2+q}{\sum\limits\_{i=p\_1+1}^{\left(1+\eta\_1^{\ell}\right)} \left(\sum\limits\_{j=1}^{q} \mathcal{C}\_{j+p-2}^{i-1} \frac{w(\left(k+j-i-1\right)}{w(\left(0\right)+2q\_{\inf}^{i(j-1)}w(\left(k\right)\right)}\right)^{2}} + \sum\limits\_{i=(n+1)/2+q+1}^{n-1} \left(\sum\limits\_{j=1}^{(n+1)/2+q-\eta} \mathcal{C}\_{p+i+j-q-2-(n+1)/2}^{i-(n+1)/2} \frac{w(\left(k+j-q-1-(n+1)/2\right)}{w(\left(0\right)+2q\_{\inf}^{i(j-1)}w(\left(k\right)\right)}\right)^{2}} \\ \mathcal{N}\_{\upsilon} = \sqrt{\sum\limits\_{i=q+1}^{(n+1)/2+q} \mathcal{C}\_{j+p-2}^{i-1} \frac{w(\left(k+j-i-1\right))}{w(\left(0\right)+2q\_{\inf}^{i(j-1)}w(\left(k\right)\right)}}^{2} + \sum\limits\_{i=(n+1)/2+q\_{\inf}+1}^{n} \left(\sum\limits\_{j=1}^{(n+1)/2+q-\eta} \mathcal{C}\_{p+i+j-q-2-(n+1)/2}^{i-(n+1)/2} \frac{w(\left(k+j-q\_{2}-1-(n+1)/2\right)}{w(\left(0\right)+2q\_{\inf}^{i(j-1)}w(\left(k\right)\right)}\right)^{2}} \end{array} \tag{16}$$

#### **Proof of Theorem 5.** Please see Appendix A.1.

**Theorem 6.** *If a raw data disturbance occurs in the second half of the original sequence, excluding the last data, then the difference between β* 1*<sup>v</sup> and β<sup>v</sup> satisfies expression (12). Let ϕ*<sup>3</sup> *be the largest integer less than* (*n* − *k*)/2*, ϕ*<sup>4</sup> *be the largest integer less than* (*n* − *k* + 1)/2*, then Mv and Nv are: (i). If sequence size n is even, then k* = *n*/2 + 1, *n*/2 + 2, . . . , *n* − 1 *and*

$$\begin{array}{l} \mathcal{M}\_{\upsilon} = \sqrt{\frac{n-q\_{4}}{i\upsilon+n/2+q\_{4}} \left( \sum\_{j=1}^{i-n/2+q\_{4}} \mathcal{C}\_{j+p-2}^{i-1} \frac{w(|k+j-i-1|)}{w(0)+2\sum\_{h=1}^{u(\cdot)-j+1} w(h)} \right)^{2} + \sum\_{i=n+1-q\_{4}}^{n-1} \left( \sum\_{j=1}^{n/2} \mathcal{C}\_{p+i+j+q-2-n}^{i+j+q-1-n} \frac{w(|k+j+q-1-n|)}{w(0)+2\sum\_{h=1}^{u(\cdot)-j} w(h)} \right)^{2}}{\sum\_{j=n/2+1-q\_{4}}^{n-q\_{4}} \left( \sum\_{j=1}^{i-n/2+q\_{4}} \mathcal{C}\_{j+p-2}^{i-1} \frac{w(|k+j-i-1|)}{w(0)+2\sum\_{h=1}^{u(\cdot)-j+1} w(h)} \right)^{2} + \sum\_{i=n+1-q\_{4}}^{n-1} \left( \sum\_{j=1}^{n/2} \mathcal{C}\_{p+i+j+q-2-n} \frac{w(|k+j+q\_{4}-1-n|)}{w(0)+2\sum\_{h=1}^{u(\cdot)-j+1} w(h)} \right)^{2}} \end{array} \tag{17}$$

*(ii). If the sequence size n is odd, then k* = (*n* + 1)/2 + 1, (*n* + 1)/2 + 2, . . . , *n* − 1 *and*

$$\begin{array}{l} \mathcal{M}\_{\nu} = \sqrt{\frac{\sum\_{i=q+1/2-qj}^{n-q\_{4}} \left( \sum\_{j=1}^{i+q\_{3}-(n-1)/2} \mathcal{C}\_{j+p-2}^{i-1} \frac{w(|k+j-i-1|)}{w(0)+2\sum\_{k=1}^{i(k+1-j)} w(k)} \right)^{2} + \sum\_{i=n-q\_{4}+1}^{n-1} \left( \sum\_{j=1}^{(n+1)/2+q\_{3}-q\_{4}-n} \mathcal{C}\_{p+i+j+q\_{4}-2-n}^{i+j+q\_{4}-n-1-i} \frac{w(|k+j+q\_{4}-1-i|)}{w(0)+2\sum\_{k=1}^{i(k+1-j)} w(k)} \right)^{2}}{\sum\_{i=-n-q\_{4}-n}^{n-1} \left( \sum\_{j=1}^{i+q\_{3}-(n-1)/2} \mathcal{C}\_{j+p-2}^{i-1} \frac{w(|k+j-i-1|)}{w(0)+2\sum\_{k=1}^{i(k+1-j)} w(k)} \right)^{2} + \sum\_{i=n-q\_{4}+1}^{n-1} \left( \sum\_{j=1}^{i(n-1)/2} \mathcal{C}\_{p+i+j+q\_{4}-2-n}^{i+j+q\_{4}-n-1-i} \frac{w(|k+j-i-1|)}{w(0)+2\sum\_{k=1}^{i(k+1-j)} w(k)} \right)^{2}} \end{array} \tag{18}$$

#### **Proof of Theorem 6.** Please see Appendix A.2.

**Theorem 7.** *If a raw data disturbance occurs in the last data of the original sequence, that is <sup>x</sup>*1(0)(*n*) = *<sup>x</sup>*(0)(*n*) + *<sup>ε</sup>, then the difference between <sup>β</sup>* 1*<sup>v</sup> and β<sup>v</sup> satisfies expression (12); however, the values of parameters Mv and Nv are determined by the size of samples.*

*(i). If sequence size n is even, then*

$$\begin{array}{l} \mathcal{M}\_{v} = \sqrt{\frac{\sum\_{i=1}^{n/2-1} \left( \sum\_{j=1}^{i} \mathcal{C}\_{i+p-j-1}^{i-j} \frac{w(n/2-j)}{w(0) + 2\sum\_{h=1}^{a(n/2+j)} w(h)} \right)^{2}}{\sum\_{i=1}^{n/2} \left( \sum\_{j=1}^{i} \mathcal{C}\_{i+p-j-1}^{i-j} \frac{w(n/2-j)}{w(0) + 2\sum\_{h=1}^{a(n/2+j)} w(h)} \right)^{2}} . \tag{19}$$

(ii). *If sequence size n is odd, then*

$$\begin{array}{l} \mathcal{M}\_{\upsilon} = \sqrt{\frac{\sum\_{i=1}^{(n-1)/2} \left( \sum\_{j=1}^{i} \mathcal{C}\_{i+p-j-1}^{i-j} \frac{w((n+1)/2-j)}{w(0) + 2\sum\_{h=1}^{a((n+1)/2+j-1)} w(h)} \right)^{2}}{\sqrt{\sum\_{i=1}^{(n+1)/2} \left( \sum\_{j=1}^{i} \mathcal{C}\_{i+p-j-1}^{i-j} \frac{w((n+1)/2-j)}{w(0) + 2\sum\_{h=1}^{a((n+1)/2+j-1)} w(h)} \right)^{2}}} . \tag{20} \end{array} \tag{20}$$

**Proof of Theorem 7.** Please see Appendix A.3.

#### **4. Perturbation Analysis of Model 2: The** *NIGM*(*p***,1**) **Model in [23]**

The fractional-order accumulate discrete grey model in Section 3 emphasizes old data in the modeling samples. To provide more weight to the new information, Wu et al. [23] provide the following model.

**Definition 5** [23]**.** *Let X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / *be the original sequence, and <sup>X</sup>*(1−*p*) = . *<sup>x</sup>*(1−*<sup>p</sup>*)(1), *<sup>x</sup>*(1−*<sup>p</sup>*)(2),..., *<sup>x</sup>*(1−*<sup>p</sup>*)(*n*) / (0 < *p* < 1) *be the fractional (1* − *p)-order reverse accumulated generating operator sequence of X*(0) *(see Definition 2). Set mean formula*

*<sup>z</sup>*(0)(*k*) = *<sup>x</sup>*(0)(*k*) + *<sup>x</sup>*(0)(*<sup>k</sup>* <sup>−</sup> <sup>1</sup>) /2 *and formula <sup>d</sup>*(1)*x*(1−*<sup>p</sup>*)(*k*) = *<sup>x</sup>*(1−*<sup>p</sup>*)(*k*) − *<sup>x</sup>*(1−*<sup>p</sup>*)(*<sup>k</sup>* − <sup>1</sup>)*, then the new information prior grey model N IGM*(*p*, 1) *is defined by:*

$$d\_{(1)}x\_{(1-p)}(k) + \mathfrak{r}\_1 z^{(0)}(k) = \mathfrak{r}\_2 \tag{21}$$

The least squares estimation of model parameters *τ*<sup>1</sup> and *τ*<sup>2</sup> in (21) can be obtained by:

$$
\begin{bmatrix} \mathfrak{r}\_1\\ \mathfrak{r}\_2 \end{bmatrix} = \left(E^T E\right)^{-1} E^T \mathcal{U} \tag{22}
$$

where

$$E = \begin{bmatrix} -z^{(0)}(2) & 1\\ -z^{(0)}(3) & 1\\ \vdots & \vdots\\ -z^{(0)}(n) & 1 \end{bmatrix}' \quad . \quad \mathcal{U} = \begin{bmatrix} d\_{(1)}\mathfrak{x}\_{(1-p)}(2) \\ d\_{(1)}\mathfrak{x}\_{(1-p)}(3) \\ \vdots \\ d\_{(1)}\mathfrak{x}\_{(1-p)}(n) \end{bmatrix}.$$

Let *τ* = [*τ*1, *τ*2] *<sup>T</sup>*, Δ*E* and Δ*U* be the disturbance items of model parameters *E* and *U*, respectively, and *ε* be the disturbance item on raw data. Assume that the matrix norm condition *<sup>E</sup>*−<sup>1</sup> 2 Δ*E* <sup>2</sup> < 1 holds for model *NIGM*(*p*, 1). Let *μ*† = *<sup>E</sup>*−<sup>1</sup> 2 *E* , *η*† = 1− *<sup>E</sup>*−<sup>1</sup> 2 Δ*E* <sup>2</sup>, *ητ* = *U* − *Eτ*. Based on Lemma 1, the perturbation analysis results of *NIGM*(*p*, 1) are summarized in the following Theorems 8 to 10. Theorem 8 directlycomes from Wu et al. [23], while Theorems 9 to 10 are derived from the corrected parameters.

**Theorem 8** [23]**.** *If the disturbance ε occurs in the first data of the original sequence, that is <sup>x</sup>*1(0)(1) = *<sup>x</sup>*(0)(1) + *<sup>ε</sup>, then the difference between solutions <sup>τ</sup>*<sup>1</sup> *and <sup>τ</sup> is denoted by δ*(*x*(0)(1)) *and it satisfies:*

$$\|\delta(\mathfrak{x}^{(0)}(1))\| \le |\mathfrak{e}| \frac{\mu\_{\mathfrak{t}}}{\eta\_{\mathfrak{t}}} \left( \frac{||\mathfrak{r}||}{2||E||} + \frac{1}{||E||} + \frac{\mu\_{\mathfrak{t}}}{\eta\_{\mathfrak{t}}} \frac{||\eta\_{\mathfrak{t}}||}{2||E||^{2}} \right). \tag{23}$$

**Theorem 9.** *If the disturbance occurs in the non-boundary nodes, that is <sup>x</sup>*1(0)(*i*) = *<sup>x</sup>*(0)(*i*) + *<sup>ε</sup>,* (*i* = 2, 3, ... , *n* − 1)*, then the difference between the solutions of model NIGM*(*p*, 1) *with and without data disturbance is denoted by δ*(*x*(0)(*i*))*, and it satisfies:*

$$\|\delta(\mathbf{x}^{(0)}(i))\| \le |\varepsilon| \frac{\mu\_{\mathsf{T}}}{\eta\_{\mathsf{T}}} \left( \frac{\sqrt{2} \|\mathsf{r}\|}{2 \|E\|} + \frac{\sqrt{\sum\_{j=2}^{i} \left(\mathbb{C}\_{j-2-p}^{j-1}\right)^{2} + 1}}{\|E\|} + \frac{\mu\_{\mathsf{T}}}{\eta\_{\mathsf{T}}} \frac{\sqrt{2} \|\eta\_{\mathsf{T}}\|}{2 \|E\|^{2}} \right). \tag{24}$$

**Proof of Theorem 9.** (i). When *<sup>x</sup>*1(0)(2) = *<sup>x</sup>*(0)(2) + *<sup>ε</sup>*, the disturbance items of parameters Δ*E* and Δ*U* are:

$$
\Delta E = \begin{bmatrix} -\varepsilon/2 & 0 \\ -\varepsilon/2 & 0 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{bmatrix}, \qquad \Delta \mathcal{U} = \begin{bmatrix} p\varepsilon \\ -\varepsilon \\ 0 \\ \vdots \\ 0 \end{bmatrix}.
$$

$$
\text{Then, } \left| \left| \Delta E \right| \right| = \sqrt{2} |\varepsilon|/2, \left| \left| \Delta \mathcal{U} \right| \right| = |\varepsilon|\sqrt{p^2 + 1}.
$$

(ii). When *<sup>x</sup>*1(0)(3) = *<sup>x</sup>*(0)(3) + *<sup>ε</sup>*, the disturbance items of model parameters <sup>Δ</sup>*<sup>E</sup>* and Δ*U* are:

$$
\Delta E = \begin{bmatrix} 0 & 0 \\ -\varepsilon/2 & 0 \\ -\varepsilon/2 & 0 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{bmatrix}, \qquad \Delta U = \begin{bmatrix} -\mathbb{C}\_{1-p}^2 \varepsilon \\ p\varepsilon \\ -\varepsilon \\ 0 \\ \vdots \\ 0 \end{bmatrix}.
$$

Then, Δ*E* <sup>2</sup> <sup>=</sup> <sup>√</sup>2|*ε*|/2, Δ*U* <sup>=</sup> <sup>|</sup>*ε*<sup>|</sup> %(*p*<sup>4</sup> <sup>−</sup> <sup>2</sup>*p*<sup>3</sup> <sup>+</sup> <sup>5</sup>*p*2)/4 <sup>+</sup> 1. Likewise, (iii) when *<sup>x</sup>*1(0)(*i*) = *<sup>x</sup>*(0)(*i*) + *<sup>ε</sup>*, (*<sup>i</sup>* <sup>=</sup> 4, ... *<sup>n</sup>* <sup>−</sup> <sup>1</sup>), the disturbance items of

model parameters Δ*E* and Δ*U* are:

$$
\Delta E = \begin{bmatrix} 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ -\varepsilon/2 & 0 \\ -\varepsilon/2 & 0 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{bmatrix}, \qquad \Delta U = \begin{bmatrix} -\mathbf{C}\_{i-p}^{i-1} - \mathbf{g}^{\varepsilon} \\ -\mathbf{C}\_{i-p}^{i-2} - \mathbf{c}^{\varepsilon} \\ \vdots \\ p\varepsilon \\ -\varepsilon \\ 0 \\ \vdots \\ 0 \end{bmatrix}.$$

Then,  $||\Delta E||\_2 = \sqrt{2}|\varepsilon|/2$ ,  $||\Delta l|| = |\varepsilon|\sqrt{\sum\_{j=2}^{i} \left(\mathcal{C}\_{j-2-p}^{j-1}\right)^2 + 1}$ .

Summizing the above three scenarios, the perturbation result (24) is obtained.  $\Box$ 

**Theorem 10.** *If the disturbance occurs in the last sample point, that is <sup>x</sup>*1(0)(*n*) = *<sup>x</sup>*(0)(*n*) + *<sup>ε</sup>, then the difference between the solutions of model NIGM*(*p*, 1) *is denoted by δ*(*x*(0)(*n*))*, and it satisfies:*

$$\|\delta(\mathbf{x}^{(0)}(n))\| \le |\varepsilon| \frac{\mu\_{\mathbf{t}}}{\eta\_{\mathbf{t}}} \left( \frac{||\mathbf{r}||}{2||E||} + \frac{\sqrt{\sum\_{j=2}^{n} \left(C\_{j-2-p}^{j-1}\right)^{2}}}{||E||} + \frac{\mu\_{\mathbf{t}}}{\eta\_{\mathbf{t}}} \frac{||\eta\_{\mathbf{t}}||}{2||E||^{2}} \right). \tag{25}$$

**Proof of Theorem 10.** When *<sup>x</sup>*1(0)(*n*) = *<sup>x</sup>*(0)(*n*) + *<sup>ε</sup>*, the disturbance items of parameters <sup>Δ</sup>*<sup>E</sup>* and Δ*U* are:

$$
\Delta E = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ -\varepsilon/2 & 0 \end{bmatrix}, \qquad \Delta U = \begin{bmatrix} -\mathcal{C}\_{n-p-2}^{n-1} \\ -\mathcal{C}\_{n-p-3}^{n-2} \\ \vdots \\ -\mathcal{C}\_{1-p}^{2} \varepsilon \\ p\varepsilon \end{bmatrix}.
$$

Then, ||Δ*E*||<sup>2</sup> = |*ε*|/2, ||Δ*U*|| = |*ε*| & ∑*n <sup>j</sup>*=<sup>2</sup> (*Cj*−<sup>1</sup> *<sup>j</sup>*−2−*p*) 2 . Apply these norms to Lemma 1; perturbation bound (25) is proved.

Now, we consider the effect of operator *Dv* on the stability of model *NIGM*(*p*, 1). *X*(0) = . *x*(0)(1), *x*(0)(2),..., *x*(0)(*n*) / and *Xv* = {*xv*(1), *xv*(2),..., *xv*(*n*)} are the same as those defined in Section 3. Set *x* (1−*p*) *<sup>v</sup>* (*k*) = <sup>∑</sup>*<sup>n</sup> <sup>i</sup>*=*<sup>k</sup> <sup>C</sup>i*−*<sup>k</sup> <sup>i</sup>*−*k*−*pxv*(*i*), *<sup>d</sup>*(1)*<sup>x</sup>* (1−*p*) *<sup>v</sup>* (*k*) = *<sup>x</sup>* (1−*p*) *<sup>v</sup>* (*k*) <sup>−</sup>

*x* (1−*p*) *<sup>v</sup>* (*<sup>k</sup>* <sup>−</sup> <sup>1</sup>) and *<sup>z</sup>* (0) *<sup>v</sup>* (*k*) = (*xv*(*k*) + *xv*(*k* − 1))/2. Let *τ<sup>v</sup>* be the least squares solution of model *NIGM*(*p*, 1) based on sequence *Xv*. Then we have

$$\pi\_{\upsilon} = \left(E\_{\upsilon}^{T} E\_{\upsilon}\right)^{-1} E\_{\upsilon}^{T} \mathcal{U}\_{\upsilon} \tag{26}$$

where

$$E\_{\upsilon} = \begin{bmatrix} -z\_{\upsilon}^{(0)}(2) & 1 \\ -z\_{\upsilon}^{(0)}(3) & 1 \\ \vdots & \vdots \\ -z\_{\upsilon}^{(0)}(n) & 1 \end{bmatrix}, \qquad \mathcal{U}\_{\upsilon} = \begin{bmatrix} d\_{(1)}x\_{\upsilon}^{(1-p)}(2) \\ d\_{(1)}x\_{\upsilon}^{(1-p)}(3) \\ \vdots \\ d\_{(1)}x\_{\upsilon}^{(1-p)}(n) \end{bmatrix}.$$

Let <sup>Δ</sup>*Ev* and <sup>Δ</sup>*Uv* be the disturbance items of *Ev* and *Uv*, respectively, and *<sup>τ</sup>*1*<sup>v</sup>* be the new model solution. Assume that *<sup>E</sup>*−<sup>1</sup> *<sup>v</sup>* 2 Δ*Ev* <sup>2</sup> < 1 holds and let *ηv*† = 1− *<sup>E</sup>*−<sup>1</sup> *<sup>v</sup>* 2 Δ*Ev* 2, *μv*† = *<sup>E</sup>*−<sup>1</sup> *<sup>v</sup>* 2 *Ev* , *ηv<sup>τ</sup>* = *Uv* − *Evτv*. The perturbation analysis results of model *NIGM*(*p*, 1) with operator *Dv* are stated in the following theorems.

**Theorem 11.** *If a raw data disturbance occurs in the first data of the original sequence, that is <sup>x</sup>*1(0)(1) = *<sup>x</sup>*(0)(1) + *<sup>ε</sup>, then the difference between <sup>τ</sup>*1*<sup>v</sup> and <sup>τ</sup><sup>v</sup> satisfies:*

$$\|\delta(\mathbf{x}^{(0)}(1))\| \le |\varepsilon| \frac{\mu\_{\upsilon\mathsf{t}}}{\eta\_{\upsilon\mathsf{t}}} \left( \frac{P\_{\upsilon} ||\mathsf{r}\_{\upsilon}||}{||E\_{\upsilon}||} + \frac{Q\_{\upsilon}}{||E\_{\upsilon}||} + \frac{\mu\_{\upsilon\mathsf{t}}}{\eta\_{\upsilon\mathsf{t}}} \frac{P\_{\upsilon} ||\eta\_{\upsilon\mathsf{r}}||}{||E\_{\upsilon}||^{2}} \right) \tag{27}$$

*where (i). If sequence size n is even, then:*

$$\begin{array}{l} P\_{\upsilon} = \frac{1}{2} \sqrt{\sum\_{i=1}^{n/2-1} \left( \frac{w(i-1)}{w(0) + 2\sum\_{k=1}^{\upsilon(i)} w(k)} + \frac{w(i)}{w(0) + 2\sum\_{k=1}^{\upsilon(i+1)} w(k)} \right)^2 + \left( \frac{w(n/2-1)}{w(0) + 2\sum\_{k=1}^{\upsilon(n/2)} w(k)} \right)^2}{n/2 - 1} \\ Q\_{\upsilon} = \sqrt{\sum\_{i=1}^{n/2-1} \left( \frac{w(i-1)}{w(0) + 2\sum\_{k=1}^{\upsilon(i)} w(k)} + \sum\_{j=1}^{n/2-1} \mathbb{C}\_{j-p-1}^{j} \frac{w(i+j-1)}{w(0) + 2\sum\_{k=1}^{\upsilon(i+j)} w(k)} \right)^2 + \left( \frac{w(n/2-1)}{w(0) + 2\sum\_{k=1}^{\upsilon(n/2)} w(k)} \right)^2} \end{array} \tag{28}$$

*(ii). If sequence size n is odd, then:*

$$\begin{array}{l} P\_{\upsilon} = \frac{1}{2} \sqrt{\sum\_{i=1}^{(n-1)/2} \left( \frac{w(i-1)}{w(0) + 2\sum\_{k=1}^{a(i)} w(k)} + \frac{w(i)}{w(0) + 2\sum\_{k=1}^{(i+1)} w(k)} \right)^2 + \left( \frac{w((n-1)/2)}{w(0) + 2\sum\_{k=1}^{a((n+1)/2)} w(k)} \right)^2}{\sum\_{i=1}^{(n-1)/2} \left( \frac{w(i-1)}{w(0) + 2\sum\_{k=1}^{a(i)} w(k)} + \sum\_{j=1}^{(n+1)/2-i} C\_{j-p-1} \frac{w(i+j-1)}{w(0) + 2\sum\_{k=1}^{a(i)+j} w(k)} \right)^2 + \left( \frac{w((n-1)/2)}{w(0) + 2\sum\_{k=1}^{a((n+1)/2)} w(k)} \right)^2} \tag{29}$$

**Proof of Theorem 11.** (i). The sequence size *<sup>n</sup>* is even. When *<sup>x</sup>*1(0)(1) = *<sup>x</sup>*(0)(1) + *<sup>ε</sup>*, according to the logic of *NIGM*(*p*, 1), we have

$$
\Delta E\_{\mathcal{V}} = \begin{bmatrix} -\frac{\xi}{2} \left( 1 + \frac{w(1)}{w(0) + 2w(1)} \right) & 0 \\ -\frac{\xi}{2} \left( \frac{w(1)}{w(0) + 2w(1)} + \frac{w(2)}{w(0) + 2\sum\_{i} w(i)} \right) & 0 \\ -\frac{\xi}{2} \left( \frac{w(n/2 - 1)}{w(0) + 2\sum\_{i=1}^{2} w(0)} + \frac{w(n/2 - 1)}{w(n/2)} \right) & 0 \\ \vdots & \vdots \\ -\frac{\xi}{2} \left( \frac{w(n/2 - 1)}{w(0) + 2\sum\_{i=1}^{2} w(0) + w(0) + 2\sum\_{i=1}^{2} w(i)} \right) & 0 \\ -\frac{\xi}{2} \underbrace{\frac{w(n/2 - 1)}{w(0) + 2w(1)}}\_{\text{ww}(0) + 2\sum\_{i=1}^{2} w(0)} & 0 \\ 0 & \frac{\xi}{2} \left( \frac{w(n/2 - 1)}{w(0) + 2w(1)} \right) & \frac{\xi}{2} \left( \frac{w(n/2 - 1)}{w(0) + 2\sum\_{i=1}^{2} w(0)} \right) \\ 0 & -\xi \frac{\frac{w(n/2 - 1)}{w(0) + 2\sum\_{i=1}^{2} w(0)}}{w(0) + 2\sum\_{i=1}^{2} w(0)} & -\xi \frac{\frac{w(n/2 - 1)}{w(0) + 2\sum\_{i=1}^{2} w(0)}}{w(0) + 2\sum\_{i=1}^{2} w(0)} \\ 0 & 0 & 0 \end{bmatrix}
$$

Then, calculate the matrix norm ||Δ*Ev*||<sup>2</sup> and vector norm ||Δ*Uv*||, extract parameter *ε* and apply Lemma 1; then, the parameter expressions in Equation (28) are obtained. (ii). When sequence size *n* is odd, we obtain:

Δ*Ev* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ − *ε* 2 1 + *<sup>w</sup>*(1) *w*(0)+2*w*(1) 0 − *ε* 2 ⎛ ⎜⎝ *<sup>w</sup>*(1) *<sup>w</sup>*(0)+2*w*(1) <sup>+</sup> *<sup>w</sup>*(2) *w*(0)+2 *α*(3) ∑ *h*=1 *w*(*h*) ⎞ ⎟⎠ 0 . . . . . . − *ε* 2 ⎛ ⎜⎝ *<sup>w</sup>*((*n*−3)/2) *w*(0)+2 *α*((*n*−1)/2) ∑ *h*=1 *w*(*h*) + *<sup>w</sup>*((*n*−1)/2) *w*(0)+2 *α*((*n*+1)/2) ∑ *h*=1 *w*(*h*) ⎞ ⎟⎠ 0 − *ε* 2 *w*((*n*−1)/2) *w*(0)+2 *α*((*n*+1)/2) ∑ *h*=1 *w*(*h*) 0 0 0 . . . . . . 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *Uv* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −*ε* − *ε* (*n*−1)/2 ∑ *i*=1 *Ci i*−*p*−1 *w*(*i*) *w*(0)+2 *α*(*i*+1) ∑ *h*=1 *w*(*h*) <sup>−</sup>*<sup>ε</sup> <sup>w</sup>*(1) *<sup>w</sup>*(0)+2*w*(1) − *ε* (*n*−3)/2 ∑ *i*=1 *Ci i*−*p*−1 *w*(*i*+1) *w*(0)+2 *α*(*i*+2) ∑ *h*=1 *w*(*h*) . . . <sup>−</sup>*<sup>ε</sup> <sup>w</sup>*((*n*−3)/2) *w*(0)+2 *α*((*n*−1)/2) ∑ *h*=1 *w*(*h*) + *εp <sup>w</sup>*((*n*−1)/2) *w*(0)+2 *α*((*n*+1)/2) ∑ *h*=1 *w*(*h*) <sup>−</sup>*<sup>ε</sup> <sup>w</sup>*((*n*−1)/2) *w*(0)+2 *α*((*n*+1)/2) ∑ *h*=1 *w*(*h*) 0 . . . 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

Then, calculate the matrix norm ||Δ*Ev*||<sup>2</sup> and vector norm ||Δ*Uv*||, extract parameter *ε* and apply Lemma 1; then, the parameter expressions in Equation (29) are obtained.

**Theorem 12.** *If a raw data disturbance occurs in the non-boundary nodes of the first half of the original sequence, that is <sup>x</sup>*1(0)(*k*) = *<sup>x</sup>*(0)(*k*) + *<sup>ε</sup>*, *then the difference between <sup>τ</sup>*1*<sup>v</sup> and <sup>τ</sup><sup>v</sup> satisfies expression (27); however, the values of parameters Pv and Qv are determined by the following scenarios. Let parameters ϕ*<sup>1</sup> *and ϕ*<sup>2</sup> *be the same as those defined in Theorem 5, then Pv and Qv are: (i). If sequence size n is even, then k* = 2, 3, . . . , *n*/2 *and*

$$\begin{split} P\_{\mathcal{V}} &= \frac{1}{2} \sqrt{\frac{\frac{w(|k-q\_{1}-1|)}{w(0)+2\sum\_{h=1}^{[0,\varphi]}w(h)}}{w(0)+2\sum\_{h=1}^{[0,\varphi]}w(h)}^{2} + \sum\_{i=q\_{1}+1}^{[1,\varphi]} \left(\frac{w(|k-i-1|)}{w(0)+2\sum\_{h=1}^{[i,i]}w(h)} + \frac{w(|k-i|)}{w(0)+2\sum\_{h=1}^{[i,i]}w(h)}\right)^{2} + \left(\frac{w(|k-n/2-q\_{1}|)}{w(0)+2\sum\_{h=1}^{[n/2,\varphi]}w(h)}\right)^{2}}{Q\_{\mathcal{V}} &= \sqrt{\sum\_{i=1}^{[1,\varphi]}\left(\sum\_{h=1}^{[i,i]}C\_{[i-p,2]}^{[i,i]}\frac{w(|k-i|)}{w(0)+2\sum\_{h=1}^{[i,i]}w(h)} + \frac{w(|k-i|)}{w(0)+2\sum\_{h=1}^{[n/2,\varphi]}w(h)}\right)^{2}} \\ &\quad \underbrace{\left(\sum\_{i=1}^{[1,\varphi]}\left(\sum\_{h=1}^{[i,i]}C\_{[i-p,2]}^{[i,i]}\frac{w(h)}{w(0)+2\sum\_{h=1}^{[i,i]}w(h)} + \frac{w(|k-i|)}{w(0)+2\sum\_{h=1}^{[i,i]}w(h)}\right)^{2}}\_{2} \right)^{2}}\_{2} \end{split}$$

$$\text{(ii). If sequence size } n \text{ is odd, then } k = 2, 3, \dots, (n+1)/2 \text{ and}$$

*Pv* = <sup>1</sup> 2 2334 *<sup>w</sup>*(|*k*−*ϕ*1−1|) *<sup>w</sup>*(0)+2∑*α*(*ϕ*1+1) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) 2 + (*n*+1)/2+*ϕ*2−1 ∑ *i*=*ϕ*1+1 *<sup>w</sup>*(|*k*−*i*−1|) *w*(0)+2∑*α*(*i*+1) *<sup>h</sup>*=<sup>1</sup> *<sup>w</sup>*(*h*) <sup>+</sup> *<sup>w</sup>*(|*k*−*i*|) *w*(0)+2∑*α*(*i*) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) 2 + *<sup>w</sup>*(|*k*−(*n*+1)/2−*ϕ*2|) *<sup>w</sup>*(0)+2∑*α*((*n*+1)/2+*ϕ*2) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) 2 , *Qv* = 2334 *ϕ*1 ∑ *i*=1 (*n*+1)/2+*ϕ*<sup>2</sup> ∑ *j*=*ϕ*1+1 *Cj*−<sup>1</sup> *j*−*p*−2 *w*(|*k*−*j*|) *w*(0)+2∑*α*(*j*) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) 2 + (*n*+1)/2+*ϕ*2−1 ∑ *i*=*ϕ*1+1 (*n*+1)/2+*ϕ*<sup>2</sup> ∑ *j*=*i*+1 *Cj*−<sup>1</sup> *j*−*p*−2 *w*(|*k*−*j*|) *w*(0)+2∑*α*(*j*) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) + *<sup>w</sup>*(|*k*−*i*|) *w*(0)+2∑*α*(*i*) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) 2 + *<sup>w</sup>*(|*k*−(*n*+1)/2−*ϕ*2|) *<sup>w</sup>*(0)+2∑*α*((*n*+1)/2+*ϕ*2) *<sup>h</sup>*=<sup>1</sup> *w*(*h*) 2 . (31)

**Proof of Theorem 12.** The proof is similar to that of Theorem 5, and we omit the details here.

**Theorem 13.** *If a raw data disturbance occurs in the second half of the original sequence, excluding the last data, let <sup>ϕ</sup>*<sup>3</sup> *and <sup>ϕ</sup>*<sup>4</sup> *be the same as those defined in Theorem 6, then the difference between <sup>τ</sup>*1*<sup>v</sup> and τ<sup>v</sup> satisfies expression (27), where parameters Pv and Qv are determined by: (i). If sequence size n is even, then k* = *n*/2 + 1, *n*/2 + 2, . . . , *n* − 1 *and*

$$\begin{split} P\_{\mathcal{D}} &= \frac{1}{2} \sqrt{\frac{w(|k+q\_{\mathcal{V}}-w/2-1|)}{w(0)+2\sum\_{i=1}^{n/2-1}w(0)}} \Big( \Big)^{2} + \sum\_{i=n/2+1}^{n-q\_{\mathcal{V}}-1} \Big( \frac{w(|k-i-1|)}{w(0)+2\sum\_{i=1}^{n(i)}w(0)} + \frac{w(|k-i|)}{w(0)+2\sum\_{i=1}^{n(i)}w(0)} \Big)^{2} + \left( \frac{w(|k+q\_{\mathcal{V}}-q\_{\mathcal{V}}|}{w(0)+2\sum\_{i=1}^{n(i)}w(0)} \right)^{2} \Big) \\ Q\_{\mathcal{V}} &= \sqrt{\sum\_{i=1}^{n/2-q\_{\mathcal{V}}} \left( \sum\_{j=n/2+1-q\_{\mathcal{S}}}^{n-q\_{\mathcal{V}}} \mathbb{C}^{j-1}\_{j-p-2} \frac{w(|k-j|)}{w(0)+2\sum\_{i=1}^{n(i)}w(h)} \right)^{2} + \sum\_{i=n/2+1}^{n-q\_{\mathcal{V}}-1} \Big( \sum\_{j=n/2}^{n-q\_{\mathcal{V}}} \mathbb{C}^{j-1}\_{j-p-2} \frac{w(|k-j|)}{w(0)+2\sum\_{i=1}^{n(i)}w(h)} + \frac{w(|k-i|)}{w(0)+2\sum\_{i=1}^{n(i)}w(h)} \Big)^{2} \\ & \qquad + \left( \frac{w(|k+q\_{\mathcal{V}}-q\_{\mathcal{V}}|}{w(0)+2\sum\_{i=1}^{n(i)}w(h)} \right)^{2} .\end{split} \tag{32}$$

*(ii). If sequence size n is odd, then k* = (*n* + 1)/2 + 1, (*n* + 1)/2 + 2, . . . , *n* − 1 *and*

$$P\_v = \frac{1}{2} \left\langle \left( \frac{w\{k + \phi\_2 - (n+1) / 2\}}{w(0) + 2\sum\_{l=1}^{w((n+1)/2 - \eta\_l)} w(l)} \right)^2 + \sum\_{l=(n+1)/2 - \eta\_l}^{n - \eta\_l - 1} \left( \frac{w([k-l-1])}{w(0) + 2\sum\_{l=1}^{w(l+1)} w(l)} + \frac{w([k-l])}{w(0) + 2\sum\_{l=1}^{w(l)} w(l)} \right)^2 + \left( \frac{w([k+\phi\_2 - (n+1) / 2])}{w(0) + 2\sum\_{l=1}^{w(n+1)} w(l)} + \frac{w([k-l-1])}{w(0) + 2\sum\_{l=1}^{w(l)} w(l)} \right)^2 \right\rangle$$

$$\begin{split} Q\_{v} &= \sqrt{\sum\_{i=1}^{s-1/2-\eta} \left( \sum\_{j=(u+1)/2-\eta\_{i}}^{s-\eta} C\_{j-\rho-2}^{\prime-1} \frac{w(\left\lVert k-j \right\rVert)}{w(0)+2\sum\_{i=1}^{u(\cdot)}w(h)} \right)^{2} + \sum\_{i=(u+1)/2-\eta\_{i}}^{s-\eta\_{i}-1} \left( \sum\_{j=i+1}^{s-\eta} C\_{j-\rho-2}^{\prime-1} \frac{w(\left\lVert k-j \right\rVert)}{w(0)+2\sum\_{i=1}^{u(\cdot)}w(h)} + \frac{w(\left\lVert k-i \right\rVert)}{w(0)+2\sum\_{i=1}^{u(\cdot)}w(h)} \right)^{2}} \\ &+ \left( \frac{w(\left\lVert k+\rho\_{u}-\rho\right\rVert)}{w(0)+2\sum\_{i=1}^{u(\cdot)-\eta\_{i}}w(h)} \right)^{2} .\end{split} \tag{33}$$

**Proof of Theorem 13.** The proof is similar to that of Theorem 6, and we omit the details here.

**Theorem 14.** *If a disturbance occurs in the last data of the original sequence: <sup>x</sup>*1(0)(*n*) = *<sup>x</sup>*(0)(*n*) + *<sup>ε</sup>, then the difference between <sup>τ</sup>*1*<sup>v</sup> and <sup>τ</sup><sup>v</sup> satisfies expression (27) with the following values of parameters Pv and Qv:*

#### *(i). If sequence size n is even, then*

$$\begin{array}{l} P\_{\mathcal{D}} = \frac{1}{2} \sqrt{\left(\frac{w(n/2-1)}{w(0) + 2\sum\_{h=1}^{w(n/2-1)} w(h)}\right)^{2} + \sum\_{i=n/2+1}^{n-1} \left(\frac{w(n-i-1)}{w(0) + 2\sum\_{h=1}^{w(i+1)} w(h)} + \frac{w(n-i)}{w(0) + 2\sum\_{h=1}^{w(i)} w(h)}\right)^{2}} \\ Q\_{\mathcal{D}} = \sqrt{\sum\_{i=1}^{n/2} \left(\sum\_{j=-n/2+1}^{n} \mathbb{C}\_{j-p-2}^{i-1} \frac{w(n-j)}{w(0) + 2\sum\_{h=1}^{w(i)} w(h)}\right)^{2} + \sum\_{i=n/2+1}^{n-1} \left(\frac{w(n-i)}{w(0) + 2\sum\_{h=1}^{w(i)} w(h)} + \sum\_{j=-\tilde{i}+1}^{n} \mathbb{C}\_{j-p-2}^{i-1} \frac{w(n-j)}{w(0) + 2\sum\_{h=1}^{w(i)} w(h)}\right)^{2}} \end{array} \tag{34}$$

*(ii). If sequence sizenis odd, then*

$$\begin{split} P\_{\mathcal{D}} &= \frac{1}{2} \sqrt{\frac{\left(\frac{w((n-1)/2)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)}\right)^{2} + \sum\_{i=\lfloor w(n+1)/2 \rfloor}^{n-1} \left(\frac{w(n-i-1)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)} + \frac{w(n-i)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)}\right)^{2}}{\sqrt{\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} \left(\sum\_{j=\lfloor w(n+1)/2 \rfloor}^{n} \frac{w(n-j)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)} + \sum\_{i=\lfloor w(n+1)/2 \rfloor}^{n-1} \frac{w(n-i)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)}} \right)^{2}} \tag{35} \\ &= \sqrt{\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} \left(\sum\_{j=\lfloor w(n+1)/2 \rfloor}^{n-1} \frac{w(n-i)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)} + \sum\_{j=\lfloor w(n+1)/2 \rfloor}^{n-1} \frac{w(n-j)}{w(0) + 2\sum\_{i=1}^{\lfloor w(n+1)/2 \rfloor} w(h)} \right)^{2}} \end{split} \tag{36}$$

**Proof of Theorem 14.** The proof is similar to that of Theorem 7, and we omit the details here.

#### **5. Numerical Performance of Operator** *Dv* **in Improving Model Stability**

*5.1. Numerical Study on Model DGMp*(1, 1)

To evaluate the effect of the fractional-order bidirectional weakening buffer operator *Dv* in improving the stability of model *DGMp*(1, 1), we consider a real case presented in [22]. This example predicts theannual cargo turnover in Jiangsu, a large coastal province in China. The original data are shown in Table 1, and the forecasting values of four compared models are presented in Table 2.

**Table 1.** The originaldata of the annual cargo turnover used in [22]. (Unit: 100 million ton-km).


**Table 2.** The forecastingresults of four different models.


\* values obtained from [22].

Table 2 reveals that the forecasting results of the *DGMp*(1, 1) model with data preprocessing based on *Dv* (*v* = 0.02) are small, which means the validity of the predicted values is maintained.

Then, data noise interferences with different noise amplitudes (ranging from −10% to 10% of the original value) are studied. The effectiveness of operator *Dv* on the prediction model is evaluated from two aspects: model stability and model accuracy. The model stability is measured by the perturbation bounds of the model parameter valueswith and without noise interference. Based on the theorems in Section 3, the perturbation bounds of the prediction models are calculated and summarized in Table 3. The model accuracy is measured by the variationsin the predicted values (VP). Let *y*0(*i*) (*i* = 1, ... *n*) be the predicted value based on the original modeling data without noise interference and *y*1(*i*) (*i* = 1, ... *n*) be the predicted value with data noise. We define *VP* = (∑*<sup>n</sup> <sup>i</sup>*=1|*y*1(*i*) − *y*0(*i*)|)/*n* and <sup>Δ</sup>*VP* <sup>=</sup> *VP*(*DGMp*(1, 1)) <sup>−</sup> *VP*(*Dv* <sup>+</sup> *DGMp*(1, 1)). If <sup>Δ</sup>*VP* falls into the positive domain, it means *Dv* improves the stability of the model *DGMp*(1, 1). The results are shown in Table 4.


**Table 3.** The perturbation bounds of four models in different noise scenarios.

\* Model index:A—*DGM*1/2(1, 1); B—*Dv* + *DGM*1/2(1, 1); C—*DGM*2/3(1, 1); D—*Dv* + *DGM*2/3(1, 1).



Table 3 shows that the *DGMp*(1, 1) models with the fractional-order bidirectional weakening buffer operator *Dv* have smaller perturbation bounds in almost all the noise scenarios except for when the data noise occurred in 2007. The results indicate that the data preprocessing method based on the fractional-order bidirectional weakening buffer operator improves the stability of the *DGMp*(1, 1) models. Table 4 shows that operator *Dv* improves the accuracy of prediction performance in most noise samples, except for the noises in the last modeling time point (in 2008). The variations in the predicted values increase with increasing noise amplitudes; however, the data preprocessing based on *Dv* effectively reduced the noise effect on the prediction results. Though *Dv* fails when the noise occurs in the last modeling time point, the data interference occurring in other time points is well controlled.

#### *5.2. Numerical Study on Model N IGM*(*p*, 1)

We now consider the effect of operator *Dv*(*v* = 0.06) in improving the stability of the model *NIGM*(*p*, 1) in forecasting the annual electricity consumption in Russia [23]. The forecasting results ofthe four comparedmodels without any data noise interference are shown in Table 5. It is clear that operator *Dv* improves the short-term (2004–2005) forecasting of *NIGM*(*p*, 1). Though its long-term (2006–2007) prediction performance is worse than those from *NIGM*(*p*, 1), its prediction accuracy is still acceptable and higher than the other two models (*DGM*(1, 1) and *GM*(0.98, 1)).


**Table 5.** The forecastingresults of *NIGM*(*p*, 1) models with and without operator *Dv*.

\* values obtained from [23].

Table 6 shows the perturbation bounds of the *NIGM*(*p*, 1) models with and without operator *Dv* in different noise scenarios. It is clear that the operator *Dv* results in smaller perturbation values in all noise scenarios and thereby improves the stability of the *NIGM*(*p*, 1) model in this numerical study. When applying the indicator *VP* defined in Section 5.1, the difference between *NIGM*(*p*, 1) with and without *Dv* can be expressed by Δ*VP* = *VP*(*NIGM*(*p*, 1)) − *VP*(*Dv* + *NIGM*(*p*, 1)). Table 7 shows the Δ*VP* results in different noise amplitude scenarios. We can see that when a data disturbance occurs in the early modeling period (2000 or 2001 in this case), operator *Dv* improves the accuracy of the prediction results in all thirty noise amplitude conditions. When noise occurs in the later modeling time point in 2002 (or 2003), the prediction without *Dv* outperforms that of with *Dv* in only four (or five) out of thirty noise perturbations. These confirm the effectiveness of operator *Dv* in improving the stability of the prediction model *NIGM*(*p*, 1).

**Table 6.** The perturbation bounds of two models in different noise scenarios.


\* Model index:E—*NIGM*(*p*, 1); F—*Dv* + *NIGM*(*p*, 1).


**Table 7.** The Δ*VP* values of *NIGM*(*p*, 1) in different noise scenarios.

The numerical cases in this section demonstrate that the fractional-order bidirectional weakening buffer operator *Dv* reduces the one-sided reaction to sample disturbances, and on the other hand, it can effectively improve the stability of time series prediction models by taking into account both the old and new information.

#### **6. Conclusions**

The performance of sequence operators in the stability of operator-based time series models has a direct impact on the validity of the model findings. Without a specific requirement for a particular part of the information in a series object, all data elements in this series should be treated equally. However, some widely used fractional-order weakening buffer operators are too subjectively biased in dealing with samples and sample disturbances. In contrast, the fractional-order bidirectional weakening buffer operator *Dv* presents a more objective approach to data processing. Both theoretical investigations and numerical experiments, based on real cases, show the favorable effects of operator *Dv*. This study reveals that the proposed fractional-order weakening buffer operator-based data preprocessing method reduces the perturbation bounds of models in data noise scenarios, thereby improving the stability and prediction accuracy of time series models. Other features of the various fractional-order operators and their effects in time series models are worth further studies.

**Author Contributions:** Conceptualization, C.L. and Y.Y.; methodology, C.L. and Y.Y.; formal analysis, C.L.; data curation, C.L. and X.Z.;writing—original draft preparation, C.L.; writing—review and editing, Y.Y.; supervision, Y.Y. and X.Z.; funding acquisition, C.L. and X.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Program of China (GrantNo. 2021YFB2601704), the Science and Technology Project of Sichuan Province (GrantNo. 2021YFS0319) and the Fundamental Research Funds for the Central Universitiesof China (GrantNo.ZHMH2022–008, J2022–062 and J2023–049).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Appendix A.1. Proof of Theorem 5*

**Proof.** (i). Sequence size *n* is even and *k* = 2, 3, ... , *n*/2. According to the proposed fractional-order operator *Dv*, the sequences with and without disturbance applied to *DGMp*(1, 1) satisfy

$$\bar{\chi}\_{\upsilon}(i) = \begin{cases} \chi\_{\upsilon}(i), & \text{when } 1 \le i \le \varphi\_1 \text{ or } n/2 + \varphi\_1 < i \le n\\ \chi\_{\upsilon}(i) + \frac{\epsilon w(|i - k|)}{\left(w(0) + 2\sum\_{j=1}^{a(i)} w(j)\right)'}, & \text{when } \varphi\_1 < i \le n/2 + \varphi\_1 \end{cases} \tag{A1}$$

Then, we can obtain

Then, calculatethe matrix norm ||Δ*Bv*||<sup>2</sup> and vector norm ||Δ*Yv*||, extract parameter *ε* and apply Lemma 1; then, the parameter expressions in Equation (15) are obtained.

(ii). Sequence size *n* is odd and *k* = 2, 3, ... ,(*n* + 1)/2. According to the fractionalorder bidirectional weakening buffer operator *Dv*, the sequences with and without disturbance applied to *DGMp*(1, 1) satisfy

$$\mathbf{x}\_{\mathbf{P}}(i) = \begin{cases} \mathbf{x}\_{\mathbf{v}}(i), \text{ when } 1 \le i \le q\_1 \text{ or } (n+1)/2 + q\_2 < i \le n\\ \mathbf{x}\_{\mathbf{v}}(i) + \frac{\text{cw}(|i-k|)}{\left(w(0) + 2\sum\_{j=1}^{a(i)} w(j)\right)}, \text{ when } q\_1 < i \le (n+1)/2 + q\_2 \end{cases} \tag{A2}$$

Then, we can obtain

Then, calculate the matrix norm ||Δ*Bv*||<sup>2</sup> and vector norm ||Δ*Yv*||, extract parameter *ε* and apply Lemma 1; then the parameter expressions in Equation (16) are obtained. The proof of Theorem 5 is now complete.

,

#### *Appendix A.2. Proof of Theorem 6*

**Proof.** (i). Sequence size *n* is even and *k* = *n*/2 + 1, *n*/2 + 2, ... , *n* − 1. According to the fractional-order bidirectional weakening buffer operator *Dv*, the sequences with and without disturbance applied to *DGMp*(1, 1) satisfy

$$\mathbf{x}\_{\upsilon}(i) = \begin{cases} \mathbf{x}\_{\upsilon}(i), \text{ when } 1 \le i < n/2 + 1 - \varrho\_4 \text{ or } n + 1 - \varrho\_4 \le i \le n\\\mathbf{x}\_{\upsilon}(i) + \frac{\varepsilon w(|i - k|)}{\left(w(0) + 2\sum\_{j=1}^{q(i)} w(j)\right)'}, \text{ when } n/2 + 1 - \varrho\_4 \le i < n + 1 - \varrho\_4 \end{cases} \tag{A3}$$

Then, we can obtain

Then, calculate the matrix norm ||Δ*Bv*||<sup>2</sup> and vector norm ||Δ*Yv*||, extract parameter *ε* and apply Lemma 1; then, the parameter expressions in Equation (17) are obtained.

(ii). Sequence size *n* is odd and *k* = (*n* + 1)/2 + 1, (*n* + 1)/2 + 2, ... , *n* − 1 . Based on operator *Dv*, the sequences with and without disturbance applied to *DGMp*(1, 1) satisfy

$$\overline{\mathbf{x}}\_{\mathbf{v}}(i) = \begin{cases} \mathbf{x}\_{\mathbf{v}}(i), & \text{when } 1 \le i < (n+1)/2 - \varrho\_3 \text{ or } n+1 - \varrho\_4 \le i \le n\\\mathbf{x}\_{\mathbf{v}}(i) + \frac{cw(|i-k|)}{\left(w(0) + 2\sum\_{j=1}^{s(i)} w(j)\right)}, & \text{when } (n+1)/2 - \varrho\_3 \le i < n+1 - \varrho\_4 \end{cases} \tag{A4}$$

Then, we can obtain

,

*Appendix A.3. Proof of Theorem 7*

Then, calculate the matrix norm ||Δ*Bv*||<sup>2</sup> and vector norm ||Δ*Yv*||, extract parameter *ε* and apply Lemma 1; then, the parameter expressions in Equation (19) are obtained. (ii). When sequence size *n* is odd, we can obtain:

Then, calculate the matrix norm ||Δ*Bv*||<sup>2</sup> and vector norm ||Δ*Yv*||, extract parameter *ε* and apply Lemma 1; then, the parameter expressions in Equation (20) are obtained. The proof of Theorem 7 is now complete.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
