A new spectral conjugate gradient method and ARIMA combined forecasting model and application

Abstract

It is of great practical significance to fit and predict actual time series. Based on the theories of time series analysis and unconstrained optimization, a new spectral conjugate gradient method–autoregressive integrated moving average combined model (FHS spectral CG–ARIMA combined model) is proposed to fit and predict the actual time series. First, combining the characteristics and advantages of different CG methods, we propose Fang–Hestenes–Stiefel algorithm (FHS). FHS satisfies the automatic descent property and has global convergence under the reasonable assumptions and Wolfe search. Second, many numerical results have been given there: compared with other related algorithms, FHS algorithm has obvious advantages. Third, FHS spectral CG–ARIMA combined model is given in detail. Fourth, the combined model is applied to fit the actual time series and the fitting effect is found to be remarkable.

Keywords

Time series analysis conjugate gradient methods FHS spectral CG-ARIMA combined model FHS algorithm global convergence

Introduction

Time series analysis is a branch of probability statistics, which has a strong applicability in signal processing, automation, information, management, financial economics, control and systems engineering, meteorological hydrology, data mining, mechanical vibration, and many other fields. From financial economics to engineering technology, from astronomy to geography and meteorology, time series analysis is used widely and has attracted much attention and is likely to be encountered in various domains.¹ In today’s rapid development of science and technology, the analysis and processing of the data has important application value, so it is more meaningful to improve the accuracy of fitting and prediction. With the rapid development of nonlinear iterative method and time series analysis, some new works have been done to realize fitting and prediction of the actual time series.

Autoregressive integrated moving average (ARIMA) model has been widely used in time series analysis. Since the ARIMA model is essentially a combination of differential operation with the autoregressive moving average model (ARMA), the parameters of the ARIMA model can be estimated by the parameters of the ARMA model.² With the extension of prediction cycle, the accuracy of conventional ARMA model is gradually reduced. Based on this, many scientific researchers have studied ways of improving the parameter estimation level of the traditional ARMA model. Shan et al.³ used two autoregressive models (AR) to estimate the parameters of ARMA model, which largely overcame the shortcomings of traditional ARMA model prediction, but its predictive accuracy is not high enough. Shan et al.^3,4 proposed the problem of estimating the parameters of ARMA model by using the method of parametric optimization estimation, but the two iterative algorithms converge slowly and their forecasting results are just passable.

For solving unconstrained optimization problems, the commonly used method is the CG method, which is especially suitable for solving large dimension problems. Its convergence rate is between Newton method and steepest descent method, and the CG method avoids the shortcoming of the Newton method to calculate the Hessen matrix, and also has a secondary termination.⁵ Spectral CG method is a variation of the classical CG method, which was first proposed by Birgin and Martinez.⁶ The main difference between it and the classical CG method is the search direction, where the spectral CG method’s search direction formula is usually $d_{0} = - g_{0}$ _, $d_{k} = - θ_{k} g_{k} + β_{k} d_{k - 1}$ , or $d_{k} = - θ_{k} g_{k} + β_{k} s_{k - 1}$ , $β_{k}$ is the conjugate parameter, and $θ_{k}$ is the spectral parameter. The spectral CG method is easy to satisfy the full descending condition or conjugate condition through the adjustment of two parameters.⁷ Numerical results show that the spectral CG algorithm is more effective than the traditional CG algorithm.⁶ So it has very important significance of theoretical research of spectral CG method. On the basis of previous studies, on the one hand, we propose a more efficient spectral CG method (FHS), which essentially promotes the combined forecasting model to improve efficiency and accuracy. On the other hand, we should pay attention to data processing as far as possible to minimize errors.

A new spectral CG method–FHS algorithm

Overview of the related methods

Consider the unconstrained optimization problem (UP)

\min f (x), x \in R^{n}

(1)

where the function

f : R^{n} \to R^{1}

is continuously differentiable.

The most commonly used method for solving this kind of problem is the CG method. Its main iterative format is

x_{k + 1} = x_{k} + α_{k} d_{k}

(2)

where

α_{k}

is the step factor, which can be determined by some methods (line search, etc.),

d_{k}

is the down search direction, and

β_{k}

is a scalar. Different CG methods are generated by different formulas of scalar parameters, and different spectral CG methods are generated according to the different search directions

d_{k}

.⁵

Some well-known CG algorithms are Fletcher–Reeves (FR) algorithm,⁸ Hestenes–Stiefel (HS) algorithm,⁹ Polak–Ribière–Polyak (PRP) algorithm,^10,11 and Dai–Yuan (DY) algorithm.¹² Among these four algorithms, FR and DY algorithm have good global convergence, while HS and PRP algorithm have fantastic numerical performance. HS algorithm for strict convex quadratic function has finite step convergence under exact line search, but for general nonstrict convex quadratic objective function, even under exact line search cannot guarantee convergent in finite steps, and global convergence cannot be guaranteed.¹³ Combining with the advantage of HS algorithm and DY algorithm, Yao et al.¹⁴ proposed an improved CG method

β_{k}^{MHS} = \frac{g_{k}^{T} (g_{k} - \frac{‖ g_{k} ‖}{‖ g_{k - 1} ‖} g_{k - 1})}{d_{k - 1}^{T} (g_{k} - g_{k - 1})}

$d_{k} = {\begin{matrix} - g_{k}, & k = 1 \\ - g_{k} + β_{k} d_{k - 1}, & k > 1 \end{matrix}$ (MHS algorithm)

And Shi et al.¹⁵ proposed a new CG method

β_{k}^{NLS} = {\begin{matrix} \frac{g_{k}^{T} (g_{k} - g_{k - 1})}{{(g_{k}^{T} d_{k - 1})}^{2} - d_{k - 1}^{T} g_{k - 1}}, & (1 - \cos θ) {‖ g_{k} ‖}^{2} > | g_{k}^{T} g_{k - 1} | \\ \max {\frac{{‖ g_{k} ‖}^{2}}{d_{k - 1}^{T} y_{k - 1}}, 0}, & otherwise \end{matrix}

$d_{k} = {\begin{matrix} - g_{k}, & k = 1 \\ - (1 + β_{k} \frac{g_{k}^{T} d_{k - 1}}{{‖ g_{k} ‖}^{2}}) g_{k} + β_{k} d_{k - 1}, & k > 1 \end{matrix}$ (NLS-DY algorithm)

The numerical performance of the MHS algorithm is between DY and HS, and each step is fully degraded and has good convergence. NLS-DY algorithm not only takes into account the numerical effect of HS and the convergence of DY to construct $β_{k}$ , but also use the modified $d_{k}$ formula proposed by Narushima et al.¹⁶

The first motivation of this paper is to combine the advantages of MHS and NLS-DY in order to provide novel algorithms with better convergence and numerical result. We choose the numerators of $β_{k}^{MHS}$ and the denominator of $β_{k}^{NLS}$ , then use Narushima’s $d_{k}$ to form our own algorithm, i.e. FHS.

FHS algorithm and its convergence analysis

For the unconstrained optimization problem (1), combining with the studies of Yao et al.¹⁴ and Shi et al.¹⁵ FHS algorithm is given here as

β_{k}^{FHS} = \frac{g_{k}^{T} (g_{k} - \frac{‖ g_{k} ‖}{‖ g_{k - 1} ‖} g_{k - 1})}{{(g_{k}^{T} d_{k - 1})}^{2} - d_{k - 1}^{T} g_{k - 1}}

(3)

d_{k} = {\begin{matrix} - g_{k}, & k = 1 \\ - (1 + β_{k} \frac{g_{k}^{T} d_{k - 1}}{{‖ g_{k} ‖}^{2}}) g_{k} + β_{k} d_{k - 1}, & k > 1 \end{matrix}

(4)

FHS algorithm implementation process:

Given an initial value $x_{1} \in R^{n}, ε > 0, d_{1} = - g_{1}, k = 1.$

Perform the Wolfe line search (WLS)

{\begin{matrix} f (x_{k} + α_{k} d_{k}) \leq f (x_{k}) + δ α_{k} g_{k}^{T} d_{k} \\ d_{k}^{T} g (x_{k} + α_{k} d_{k}) \geq σ d_{k}^{T} g_{k} \end{matrix}

(5)

where

0 < δ < σ < \frac{1}{4},

and are real numbers. From equation (5), we get

α_{k},

and according to equation (1), we obtain x_k_+1. Then calculate

f_{k + 1} = f (x_{k + 1}), g_{k + 1} = g (x_{k + 1}) .

3. If $‖ g_{k + 1} ‖ < ε$ , the minimum value is $x_{k + 1}$ ; but if $‖ g_{k + 1} ‖ \geq ε,$ go to the next step.

4. Calculate equations (3) and (4).

5. Put $k = k + 1,$ and turn to step 2.

Global convergence of FHS algorithm:

Assumptions:¹³

The level set $Ω = {x \in R^{n} | f (x) \leq f (x_{0})}$ is bounded, where $x_{0}$ is the initial point.

The function $f$ is continuously differentiable in a neighborhood $Φ$ of $Ω,$ and the function gradient satisfies the Lipschitz continuity condition, i.e. there exists a positive constant $L$ such that the following holds

‖ g (x_{1}) - g (x_{2}) ‖ < L ‖ x_{1} - x_{2} ‖, \forall x_{1}, x_{2} \in Φ

Theorem 2.1. If $g_{k} \neq 0,$ the directions generated by FHS algorithm is descendant, i.e. $\forall k \geq 1, g_{k}^{T} d_{k} < 0$ .

Proof. While $k = 1, d_{k} = - g_{k}, g_{k}^{T} d_{k} = - {‖ g_{k} ‖}^{2};$ While

k > 1, d_{k} = - (1 + β_{k} \frac{g_{k}^{T} d_{k - 1}}{{‖ g_{k} ‖}^{2}}) g_{k} + β_{k} d_{k - 1},

g_{k}^{T} d_{k} = - (1 + β_{k} \frac{g_{k}^{T} d_{k - 1}}{{‖ g_{k} ‖}^{2}}) {‖ g_{k} ‖}^{2} + β_{k} g_{k}^{T} d_{k - 1} = - {‖ g_{k} ‖}^{2},

So it is easy to verify that equation (3) satisfies the descending direction.

Lemma 2.1.¹⁷ Suppose that f satisfies the above premises, $x_{k}$ from equation (2), d_k from equation (4), $α_{k}$ satisfies the Wolfe line search (5), then Zoutendijk holds

\sum_{k = 1}^{\infty} \frac{{(g_{k}^{T} d_{k})}^{2}}{{‖ d_{k} ‖}^{2}} < + \infty

(6)

Theorem 2.2. Assuming that assumptions (1) and (2) is established, consider the conjugate gradient method based on equations (2) and (4), and $β_{k}$ from equation (3), the following holds

\lim_{k \to \infty} \inf ‖ g_{k} ‖ = 0

(7)

Proof. (Reduction to absurdity) First of all, assume that the conclusion is not established, then $\forall k > 0, \exists ε > 0$ is a real constant, and $‖ g_{k} ‖ > ε$ holds.

According to equation (4) and $g_{k}^{T} d_{k} = - {‖ g_{k} ‖}^{2}$ and let $l_{k} = 1 + β_{k} \frac{g_{k}^{T} d_{k - 1}}{{‖ g_{k} ‖}^{2}}$ , while k>0, equation (4) is squared, taken norm, and simplified as follows

\begin{array}{l} {‖ d_{k} ‖}^{2} = {(β_{k})}^{2} {‖ d_{k - 1} ‖}^{2} - 2 l_{k} d_{k}^{T} g_{k} - l_{k}^{2} {‖ g_{k} ‖}^{2} \\ = {[\frac{g_{k}^{T} (g_{k} - \frac{‖ g_{k} ‖}{‖ g_{k - 1} ‖} g_{k - 1})}{{(d_{k}^{T} d_{k - 1})}^{2} - d_{k - 1}^{T} g_{k - 1}}]}^{2} {‖ d_{k - 1} ‖}^{2} \\ - 2 l_{k} d_{k}^{T} g_{k} - l_{k}^{2} {‖ g_{k} ‖}^{2} \end{array}

then we have

\begin{array}{l} \frac{{‖ d_{k} ‖}^{2}}{{(g_{k}^{T} d_{k})}^{2}} = {[\frac{g_{k}^{T} (g_{k} - \frac{‖ g_{k} ‖}{‖ g_{k - 1} ‖} g_{k - 1})}{{(g_{k}^{T} d_{k - 1})}^{2} - d_{k - 1}^{T} g_{k - 1}}]}^{2} \frac{{‖ d_{k - 1} ‖}^{2}}{{(g_{k}^{T} d_{k})}^{2}} \\ - \frac{2 l_{k} d_{k}^{T} g_{k}}{{(g_{k}^{T} d_{k})}^{2}} - \frac{l_{k}^{2} {‖ g_{k} ‖}^{2}}{{(g_{k}^{T} d_{k})}^{2}} \\ \leq \frac{{[g_{k}^{T} (g_{k} - \frac{‖ g_{k} ‖}{‖ g_{k - 1} ‖} g_{k - 1})]}^{2}}{{(g_{k - 1}^{T} d_{k - 1})}^{2}} \frac{{‖ d_{k - 1} ‖}^{2}}{{‖ g_{k} ‖}^{4}} - \frac{{(l_{k} - 1)}^{2}}{{‖ g_{k} ‖}^{2}} \\ + \frac{1}{{‖ g_{k} ‖}^{2}} \leq \frac{{(‖ g_{k} ‖ + \frac{‖ g_{k} ‖}{‖ g_{k - 1} ‖} ‖ g_{k} ‖ ‖ g_{k - 1} ‖)}^{2}}{{(g_{k - 1}^{T} d_{k - 1})}^{2}} \frac{{‖ d_{k - 1} ‖}^{2}}{{‖ g_{k} ‖}^{4}} \\ + \frac{1}{{‖ g_{k} ‖}^{2}} = \frac{4 {‖ d_{k - 1} ‖}^{2}}{{(g_{k - 1}^{T} d_{k - 1})}^{2}} + \frac{1}{{‖ g_{k} ‖}^{2}}, \\ \frac{{‖ d_{k} ‖}^{2}}{{(g_{k}^{T} d_{k})}^{2}} \leq \sum_{i = 0}^{k} \frac{4^{k - i}}{{‖ g_{i} ‖}^{2}} \leq \sum_{i = 0}^{k} \frac{4^{k - i}}{<^{2}} = \frac{1}{<^{2}} \sum_{i = 0}^{k} 4^{k - i} \\ = \frac{1}{<^{2}} \frac{1 - 4^{k - i}}{1 - 4} = \frac{4^{k + 1} - 1}{3 <^{2}} \end{array}

It can be delivered as follows

\begin{matrix} {\frac{(g_{k}^{T} d_{k})}{{‖ d_{k} ‖}^{2}}}^{2} \geq \frac{3 ε^{2}}{4^{k + 1} - 1}, \\ \sum_{k = 1}^{\infty} \frac{{(g_{k}^{T} d_{k})}^{2}}{{‖ d_{k} ‖}^{2}} \geq \sum_{k = 1}^{\infty} \frac{3 ε^{2}}{4^{k + 1} - 1} = 3 ε^{2} \sum_{k = 1}^{\infty} \frac{1}{4^{k + 1} - 1} \end{matrix}

It is contrary to condition (6) of Lemma 2.1. If condition (6) holds, then Theorem 2.2. holds, and FHS algorithm has a global convergence.

Numerical experiments

In this section, we will use some test functions of More et al.¹⁸ to test the numerical performance of FHS, NLS, MHS, and HS algorithms. In all methods, the stepsize $α_{k}$ is yielded by WLS (5), and the parameters in WLS are chosen to be $δ = 0.01, σ = 0.1$ . All codes were written in the MATLAB 2010b, and run on the computer with Intel(R) Core(TM) i5–5200U CPU @2.2.GHz and 4.00 GB SDRAM. The numerical results of FHS, NLS, MHS, and HS algorithms are shown in Tables 1 and 2, where Problem denotes the name of the test function, n denotes the dimension of the test problem, NI/NF/NG denotes the number of iterations, function evaluations, and gradient evaluations, respectively. t denotes the computing time of CPU for computing the corresponding test problem (unit: second) and f* is the optimal function value. We stop the iteration if one of the following conditions is satisfied: (i) ${‖ g_{k} ‖}_{\infty} \leq 10^{- 6}$ ; (ii) the number of iteration $NI > 10000$ ; (iii) an uphill search direction occurs. If condition (ii) or (iii) occurs, the method is deemed to fail for solving the corresponding test problem and is denoted by “err”.

Table 1.

Iterative comparison of four algorithms: FHS, NLS, MHS, HS.

Problem	n	FHS			NLS			MHS			HS
Problem	n	NI	NF	NG	NI	NF	NG	NI	NF	NG	NI	NF	NG
Rosenbrock	2	161	585	746	197	714	911	155	560	715	26	169	195
Freudenstein	2	13	50	63	14	54	68	13	46	59	10	36	46
Brown	2	7	20	27	7	20	27	9	27	36	7	19	26
Beale	2	14	32	46	16	38	54	14	34	48	16	33	49
Jennrich1	2	7	24	31	7	25	32	5	16	21	7	63	70
Jennrich2	2	9	36	45	9	37	46	err	err	err	err	err	err
Jennrich3	2	7	25	32	7	25	32	7	25	32	7	35	42
Helical	3	54	171	225	52	149	201	48	153	201	39	111	150
Box1	3	159	469	628	70	204	274	99	288	387	53	150	203
Box2	3	13	28	40	14	32	46	13	30	43	11	26	37
Box3	3	4	9	13	4	9	13	4	9	13	6	17	23
Powell	4	40	107	147	56	154	210	56	156	212	63	172	235
Wood	4	39	127	166	35	119	154	33	116	149	err	err	err
Kowalik	4	40	114	154	14	38	52	22	60	82	23	59	82
Dennis1	4	48	156	204	44	154	198	44	153	197	42	155	197
Dennis2	4	43	145	188	44	146	190	43	147	190	56	191	247
Dennis3	4	68	230	298	92	309	401	52	171	223	59	199	258
Osborne1	5	29	88	117	24	74	98	29	89	118	err	err	err
Biggs	6	117	386	503	28	88	116	76	241	317	38	117	155
Osborne2	11	10	30	40	err	err	err	err	err	err	err	err	err

Table 2.

Run time and function value comparison of four algorithms: FHS, NLS, MHS, HS.

Problem	n	FHS		NLS		MHS		HS
Problem	n	t	f*	t	f*	t	f*	t	f*
Rosenbrock	2	16.27	4.25e-009	21.89	4.83e-009	18.23	5.12e-009	4.11	1.08e-012
Freudenstein	2	1.91	48.98	2.43	48.98	2.05	48.98	1.49	48.98
Brown	2	0.54	3.03e-015	0.58	1.75e-015	0.79	5.19e-015	0.47	2.06e-016
Beale	2	1.93	1.14e-013	2.21	2.02e-016	1.94	4.19e-018	1.82	3.03e-014
Jennrich1	2	2.72	124.89	2.97	124.83	2.30	124.81	6.73	124.45
Jennrich2	2	2.76	19.32	3.02	19.30	err	err	err	err
Jennrich3	2	0.98	0.27	0.86	0.27	0.95	0.27	1.20	0.27
Helical	3	12.74	4.93e-015	10.95	5.56e-012	10.64	4.31e-017	8.13	2.54e-011
Box1	3	46.25	1.38e-011	17.69	1.18e-014	25.26	4.93e-019	12.14	3.42e-012
Box2	3	1.62	2.45e-013	2.16	8.94e-021	1.69	3.64e-014	1.70	2.91e-015
Box3	3	0.29	0	0.43	4.93e-032	0.27	4.93e-032	0.65	1.11e-031
Powell	4	6.01	1.83e-008	9.46	1.76e-008	10.08	1.95e-008	8.67	4.07e-008
Wood	4	5.75	1.10e-016	5.82	1.90e-009	5.71	1.90e-008	err	err
Kowalik	4	15.69	7.41e-012	4.49	5.93e-012	8.13	2.24e-012	7.24	4.88e-019
Dennis1	4	35.08	1.05e-005	28.06	1.05e-005	32.62	1.05e-005	27.67	1.05e-005
Dennis2	4	24.42	3.23e-007	2.07	3.23e-007	24.97	3.23e-007	27.90	3.23e-007
Dennis3	4	27.58	9.00e-021	31.67	2.29e-020	22.28	2.57e-020	19.94	3.20e-019
Osborne1	5	13.41	6..62e-006	10.34	6.62e-006	14.50	6.63e-006	err	err
Biggs	6	77.42	1.34e-019	12.95	9.16e-014	42.33	2.50e-014	17.85	1.59e-008
Osborne2	11	13.7	0.001	err	err	err	err	err	err

From Tables 1 and 2, we can see that the NI, NF, and NG of most test functions, which are calculated by FHS algorithm, are obviously less than by NLS and MHS, and t of FHS is obviously lower than NLS and MHS. So the new iterative method (FHS) is more efficient; moreover, because the superior numerical results of the HS algorithm are well known. And calculating about half of the functions, FHS algorithm performs slightly better than HS. The reduction of the number of iterations and the running time reflects the strong convergence of the FHS algorithm, where the decrease of the error indicates that the algorithm has better numerical results.

To sum up, FHS algorithm is more effective for solving unconstrained problems. And this is the basis for making combined forecasting model more efficient and accurate.

FHS spectral CG-ARIMA combined model

ARIMA(p,d,q) model, which is the autoregressive integrated moving average model, can be understood as a combination of differential operations and ARMA(p,q) model.² In order to improve the prediction accuracy of the ARMA(p,q) model, in this section, estimation of the model parameters are transformed into a unconstrained optimization problem, and then combining with the CG method proposed in the previous section (i.e. FHS), we propose an improved spectral CG-ARIMA combined model.

The determination of the objective function

With reference to the existing literature,¹⁹ ARMA(p,q) model was transformed into unconstrained optimization problems, whose model structure is as follows

\begin{array}{l} x_{t} - ϕ_{1} x_{t - 1} - \dots - ϕ_{p} x_{t - p} = ε_{t} - θ_{1} ε_{t - 1} - \dots - θ_{q} ε_{t - q} \\ x_{t} = ϕ_{1} x_{t - 1} + \dots + ϕ_{p} x_{t - p} + ε_{t} - θ_{1} ε_{t - 1} - \dots - θ_{q} ε_{t - q} . \end{array}

The above formula could be described as $x_{t} = X_{t}^{T} ω + ε$ , where

\begin{array}{l} X_{t}^{T} = [x_{t - 1}, x_{t - 2}, \dots, x_{p}, ε_{t - 1}, ε_{t - 2}, \dots, ε_{t - q}] \\ ω = {[ϕ_{1}, ϕ_{2}, \dots, ϕ_{p}, - θ_{1}, - θ_{2}, \dots, - θ_{q}]}^{T} \\ k = m = p + q, f (X_{t}, ω) = X_{t}^{T} ω \end{array}

Because the nonstationary time series can become the stationary time series by difference operation, the ARMA model can only be discussed here. For the nonlinear relations between $ω$ and $X_{t}$ of the ARMA model, we define the objective function as the sum of squares of residuals

S (ω) = \sum_{k = p + 1}^{n} ε_{t}^{2} = \sum_{k = p + 1}^{n} {[x_{t} - f (X_{t}, ω)]}^{2}

Then the estimating value of $ω$ is converted to seeking the optimal solution $ω^{*}$ , which makes $S (ω)$ achieve the minimum.

The determination of the initial value

Here, we use the long auto-regressive (AR) model, where its calculation principle starts from the equivalent system transfer function of the model, introduces the concept of inverse function, and combines the method of undetermined coefficients, then we get

{\begin{array}{l} ϕ_{1}^{0} = θ_{1}^{0} + I_{1}, \\ ϕ_{2}^{0} = θ_{2}^{0} - θ_{1}^{0} I_{1} + I_{2}, \\ ϕ_{3}^{0} = θ_{3}^{0} - θ_{2}^{0} I_{1} - θ_{1}^{0} I_{2} + I_{3}, \\ \dots, \\ ϕ_{p}^{0} = - θ_{q}^{0} I_{p - q} - \dots - θ_{2}^{0} I_{p - 2} - θ_{1}^{0} I_{p - 1} + I_{p}, \\ 0 = - θ_{q}^{0} I_{κ - q} - \dots - θ_{2}^{0} I_{κ - 2} - θ_{1}^{0} I_{κ - 1} + I_{κ}, (κ > p) \end{array}

(8)

The above equation set (8) is regarded as the p-linear equations for $ϕ_{j}^{0}$ , where $θ_{j}^{0}$ is the coefficient. When the coefficient is known, we could obtain $ϕ_{j}^{0}$ . According to the structure of p-linear equations, we know that $θ_{j}^{0} = 0$ , when $j > p$ , and the equation set changes as follows

[\begin{matrix} ϕ_{1}^{0} \\ ϕ_{2}^{0} \\ ⋮ \\ ϕ_{p - 1}^{0} \\ ϕ_{p}^{0} \end{matrix}] = [\begin{matrix} θ_{1}^{0} \\ θ_{2}^{0} \\ ⋮ \\ θ_{p - 1}^{0} \\ θ_{p}^{0} \end{matrix}] + [\begin{matrix} 1 & 0 & 0 & \dots & 0 \\ - θ_{1}^{0} & 1 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ - θ_{p - 2}^{0} & - θ_{p - 3}^{0} & - θ_{p - 4}^{0} & \dots & 0 \\ - θ_{p - 1}^{0} & - θ_{p - 2}^{0} & - θ_{p - 3}^{0} & \dots & 1 \end{matrix}] [\begin{matrix} I_{1} \\ I_{2} \\ ⋮ \\ I_{p - 1} \\ I_{p} \end{matrix}]

(9)

when

κ

varies from

p + 1

p + q

, the above equation set (9) can be changed as

[\begin{matrix} I_{p + 1} \\ I_{p + 2} \\ ⋮ \\ I_{p_{0}} \end{matrix}] = [\begin{matrix} I_{p} & I_{p - 1} & I_{p - 2} & \dots & I_{p + 1 - q} \\ I_{p + 1} & I_{p} & I_{p - 1} & \dots & I_{p + 2 - q} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ I_{p_{0} - 1} & I_{p_{0} - 2} & I_{p_{0} - 3} & \dots & I_{p} \end{matrix}] [\begin{matrix} θ_{1}^{0} \\ θ_{2}^{0} \\ ⋮ \\ θ_{q}^{0} \end{matrix}]

(10)

We can obtain the value of $θ_{j}^{0}$ by solving the above equation set (10) and after substituting in the equation set (9), the result should be iterated to get the value of $ϕ_{j}^{0}$ . According to this method, we could achieve the original value $ω_{0}$ . The value of $ω_{0}$ varies around the extreme point from which we could acquire the result after iterating several times.

Combining with the idea of conjugate direction in nonlinear programming, an improved spectral CG method (FHS algorithm) was also proposed to optimize the parameters of the ARMA(p,q) model.

Implementation process

Spectral CG-ARIMA combined model implementation process:

Given an initial value

ω_{0} = {[ω_{1}^{0}, ω_{2}^{0}, \dots, ω_{m}^{0}]}^{T} = {[ϕ_{1}^{0}, \dots, ϕ_{p}^{0}, - θ_{1}^{0}, \dots, - θ_{q}^{0}]}^{T}

allowable error $ε > 0$ , the gradient expression of the objective function $S (ω)$ at $ω_{0}$ is

g^{0} = {[g_{1}^{0}, g_{2}^{0}, \dots, g_{m}^{0}]}^{T}

where

\begin{array}{l} g_{i}^{0} = \frac{\partial S (ω)}{\partial ω_{i}} |_{ω = ω_{0}} \\ = - 2 \sum_{t = p + 1}^{N} [x_{t} - f (X_{t}, ω)] \frac{\partial f (X_{t}, ω)}{\partial ω_{i}} |_{ω = ω_{0}} \\ = - 2 \sum_{t = p + 1}^{N} ε_{t}^{0} v_{i t}^{0} \end{array}

and

\begin{array}{l} ε_{t}^{0} = x_{t} - f (X_{t}, ω) |_{ω = ω_{0}} \\ v_{i t}^{0} = \frac{\partial f (X_{t}, ω)}{\partial ω} |_{ω = ω_{0}} (i = 1, 2, \dots, m) \end{array}

simplified as

g_{i}^{0} = - 2 {(v_{0})}^{T} ε_{0}

. Thus obtaining

S_{0} = S (ω_{0})

and

d_{1} = - g_{1}, k = 1

2. Perform the Wolfe line search

{\begin{matrix} S (ω_{k} + α_{k} d_{k}) \leq S (ω_{k}) + δ α_{k} g_{k}^{T} d_{k} \\ d_{k}^{T} g (ω_{k} + α_{k} d_{k}) \geq σ d_{k}^{T} g_{k} \end{matrix}

(11)

where

0 < δ < σ < \frac{1}{3}

, and are real numbers. From equation (11) we get

α_{k}

, and according to equation (2), we obtain

ω_{k + 1} = ω_{k} + α_{k} d_{k}

. Then

S_{k + 1} = S (ω_{k + 1}), g_{k + 1} = g (ω_{k + 1})

are calculated.

3. If $‖ g_{k + 1} ‖ < ε$ , the minimum value is $ω_{k + 1}$ ; but if $‖ g_{k + 1} ‖ \geq ε$ , go to the next step.

4. If $k = p, ω_{0} = ω_{k + 1}, S_{0} = S_{k + 1}, g_{0} = g_{k + 1}, d_{0} = - g_{0}, k = 0$ , turn to step (2); otherwise, go to the next step.

5. Calculate equations (3) and (4).

6. Put $k = k + 1$ , and if ${(d_{k + 1})}^{T} g_{k + 1} \geq 0$ , repeat the operation of step (4): $ω_{0} = ω_{k + 1}$ , $S_{0} = S_{k + 1}$ , $g_{0} = g_{k + 1}, d_{0} = - g_{0}, k = 0$ , then turn to step (2); otherwise, turn to step (2) directly.

Application of FHS spectral CG-ARIMA combined model

In order to test the practicality of this model, the time series data of Table 3 in accordance with the AR (4) model were used to estimate the unknown parameters. The initial value is calculated according to the mentioned principle:

\begin{array}{l} ω_{0} = {[ϕ_{1}^{0}, ϕ_{2}^{0}, ϕ_{3}^{0}, ϕ_{4}^{0}]}^{T} \\ = {[ω_{1}^{0}, ω_{2}^{0}, ω_{3}^{0}, ω_{4}^{0}]}^{T} = {[0.5, - 0.5, 0.5, - 0.5]}^{T} \end{array}

setting parameters:

δ = 0.01, σ =0 .1, ε {=10}^{- 6}

Table 3.

Time series.

1	2	3	4	5	6	7	8	9	10
0	0.7624	0.7753	2.1032	0.4397	−0.0109	1.6246	1.0087	1.0591	−0.4130
0.0878	−0.2615	0.6624	−0.2819	−0.0199	1.6926	0.9126	0.6219	1.7362	−2.2753
−0.1077	1.0461	0.3339	0.9674	0.6129	−0.3133	−0.2458	0.8536	2.5095	−0.1610
−0.5698	−0.6434	−1.8188	−1.3315	−0.6810	−0.3389	1.1811	1.1338	0.5663	1.3288
1.1010	0.3453	−1.1152	−0.3031	−0.3883	1.4303	0.7485	0.8044	0.4499	0.8382
−0.0697	0.8779	1.6640	0.9181	−0.7974	1.0448	0.9500	−0.1540	0.2731	0.3706
−1.1659	−0.7111	−0.6395	−0.1508	0.8148	−0.6196	−1.5186	−0.0935	0.0654	1.6739
1.1931	−0.9836	−0.5880	−0.0836	−0.5598	0.8238	2.3585	1.3246	1.7612	0.4348
0.7364	1.0715	1.1936	1.6707	1.1838	−0.1052	1.1675	−0.7770	0.7645	0.5663
−0.5328	−0.7563	−0.0932	−1.1172	0.3946	0.6220	1.2726	1.1808	0.1369	−1.0471
−0.3904	−0.0038	−0.7360	−0.7829	−1.9687

The program is written in the MATLAB 2010b, and run on the computer with Intel(R) Core(TM) i5–5200U CPU @2.2.GHz and 4.00 GB SDRAM, the results after five iterations is:

ω_{5} = {[1.0000, - 0.3934, 0.0106, - 0.0472]}^{T}

and the corresponding model structure is

x_{t} = x_{t - 1} - 0.3934 x_{t - 2} + 0.0106 x_{t - 3} - 0.0472 x_{t - 4} + ε_{t}

In Figure 1, the x-axis is the sampling point and the y-axis is the amplitude; in Figure 2, the x-axis is the retardation value and the y-axis is the normalized value. As can be seen from these figures, we get more accurate parameters from FHS spectral CG-ARIMA combined model. And from this model we can estimate the 106th and subsequent data.

Figure 1.

Raw signal.

Figure 2.

Auto-correlation function of prediction error.

Conclusion

In this paper, the combined model is one effective combination with a new spectral CG method and ARIMA model. Making full use of the advantages of MHS and NLS-DY, a new hybrid algorithm (FHS), which is more competitive, is formed. We give the detailed algorithm steps, prove that FHS is sufficient descent and global convergent under Wolfe search, and verify the high efficiency of FHS algorithm by several examples. ARIMA(p,d,q) model is a combination of differential operations and ARMA(p,q) model; FHS spectral CG-ARIMA combined model is an optimal theoretical estimation method of ARMA model parameters, the results of which show that the estimated parameters are more accurate and effective, and the combined model can be used for effective prediction.

In the future research, we will develop more effective hybrid gradient methods to solve large-scale unconstrained optimization problems; on the other hand, the accuracy of ARMA model parameters can be improved and used in the fitting and prediction of realistic problems.

Footnotes

Acknowledgements

The authors are very grateful to the anonymous referees for their valuable comments and useful suggestions, which improved the quality of this paper.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by the National Nature Science Foundation of China (Nos. 51405424, 51675461, 11673040).

References

Zhang

Lei

Feng

YQ.

Application of Matlab in time series analysis (Chinese). Shaanxi: Xi’an University of Electronic Science and Technology Press, 2007.

Wang

Applied time series analysis (Chinese). Beijing: China Renmin University Press, 2005.

Shan

Shi

Liu

Forecasting model for stock data based on new wavelet and conjugated-ARIMA combination model. J Comput Inform Syst 2014; 10: 3023–3030.

Shan

Liu

An improved conjugate gradient parameter estimation method for the differential autoregressive moving average model. Henan Univ Sci Technol 2015; 36: 85–90.

Cheng

Tang

YH.

Optimization theory and method (Chinese). Beijing: National Defense Industry Press, 2008.

Chen

Jiang

XZ.

A spectral conjugate gradient method for unconstrained optimization. J Yulin Normal Univ 2001; 43: 117–128.

Jian

Jiang

Yin

JH.

Research progress of nonlinear conjugate gradient method. J Yulin Normal Univ 2016; 37: 3–10.

Fletcher

Reeves

Function minimization by conjugate gradients. Comput J 1964; 7: 149–154.

Hestenes

Stiefel

Method of conjugate gradient for solving linear equations. J Res Natl Bureau Standard 1952; 49: 409–436.

10.

Polyak

Ribiere

Note sur la convergence de methode de directions conjuguees. Rev Franaise Informat Recherche 1969; 16: 35–63.

11.

Polyak

The conjugate gradient method in extremem problems. USSR Comp Math Math 1969; 9: 94–112.

12.

Dai

Yuan

A nonlinear conjugate gradient with a strong global convergence property. SIAM J Optim 1999; 10: 177–182.

13.

Gibert

Nocedal

Global convergence properties of conjugate gradient methods for optimization. SIAM J Optim 1992; 2: 21–42.

14.

Yao

Wei

Huang

A note about Wyl’s conjugate gradient method and its application. Appl Math Comput 2007; 191: 381–388.

15.

Shi

Shan

Liu

The conjugate gradient algorithm with perturbation factor and its convergence. J Henan Univ Sci Technol 2013; 34: 83–87.

16.

Al-Baali

Narushima

Yabe

A family of three-term conjugate gradient methods with sufficient descent property for unconstrained optimization. Comput Optim Appl 2015; 60: 89–110.

17.

Zoutendijk

Nonlinear programming computational methods. Integer Nonlin Program 1970; 143: 37–68.

18.

Garbow

Hilstrome

KE.

Testing unconstrained optimization software. ACM Trans Math Softw 1981; 7: 17–41.

19.

Zheng

Shan

An optimum algorithm for nonlinear time series ARMA parameter estimation. J Jiamusi Univ 2008; 26: 169–171.