Estimation of parameters of autoregressive models with fractional differences in the presence of additive noise

Introduction

To describe processes of various nature, equations with derivatives are increasingly used. and differences of fractional order. Despite the lack of a simple interpretation, which give derivatives, integrals and differences of integers, models described by fractional-order equations, make it possible to accurately simulate many processes in physics and technology [1$–$4]. In connection with the active development and application of equations with differences and fractional derivatives for modeling and forecasting problems, methods for estimating systems have also begun to actively develop, describing fractional-order equations and differences.

Autoregressions with fractional differences are widely used in the analysis of time series with long memory [5; 6]. There are a large number of different models with generalizations of fractional differences, such as Gegenbauer autoregressive moving average (GARMA) [7; 8], fractional ARUMA [9], seasonal autoregressive fractionally integrated moving average (SARFIMA) [10; 11], and autoregressive tempered fractionally integrated moving average (ARTFIMA) [12; 13]. Various aspects of using fractional differences for time series analysis have been considered [14; 15].

It is known that for autoregressive models of an integer order, autoregressive models with additive noise can exceed the accuracy of ARMA models and autoregressive models [16]. An overview of methods for estimating integer-order autoregressions in the presence of noise is presented in [17]. In the articles [8; 18; 19], the author considered the estimation of autoregressions with fractional-order differences in the presence of noise with a known noise ratio.

The article presents a new method for estimating parameters autoregressive models with fractional differences in the presence of additive noise with an unknown variance of additive noise.

1 Basic results

 Time series, is described by linear stochastic equations with fractional order differences:

                                                      $z_i\mathrm{<mo>=</mo>}{\displaystyle \sum _{m\mathrm{<mo>=</mo>1}}^r}{b}^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}{\Delta }^{{\alpha }_m}{z}_{i-1}+{\zeta }_i\mathrm{,}{y}_i\mathrm{<mo>=</mo>}{z}_i+{\xi }_i\mathrm{,}$                                                            (1)

 where $b^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}$ are constant coefficients; $\mathrm{0<mo><</mo>}{\alpha }_1\dots \mathrm{<mo><</mo>}{\alpha }_r\mathrm{<mo>;</mo>}$ $\Gamma \mathrm{<mo\; stretchy="false">(</mo>}\alpha \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \underset{0}{\overset{\infty }{\int }}}{e}^{-t}{t}^{\alpha -1}dt$ ;

${\Delta }^{{\alpha }_m}{z}_i\mathrm{<mo>=</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>0}}^i}{\mathrm{<mo\; stretchy="false">(</mo>}-\mathrm{1<mo\; stretchy="false">)</mo>}}^j\left(\begin{array}{c}{\alpha }_m\\ j\end{array}\right)z_{i-j}$ is fractional difference;

$\left(\begin{array}{c}{\alpha }_m\\ j\end{array}\right)\mathrm{<mo>=</mo>}\frac{\Gamma \mathrm{<mo\; stretchy="false">(</mo>}{\alpha }_m+\mathrm{1<mo\; stretchy="false">)</mo>}}{\Gamma \mathrm{<mo\; stretchy="false">(</mo>}j+\mathrm{1<mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>}{\alpha }_m-j+\mathrm{1<mo\; stretchy="false">)</mo>}}$ is generalized binomial coefficients.

It is required to estimate the unknown coefficients of the dynamic system described by (1) from the observed sequence $\left\{y_i\right\}$ with noise for the known orders $r$ , ${\alpha }_m$.

If $r$ and, ${\alpha }_m$ are unknown, it is necessary to apply algorithms based on global optimization, such as genetic algorithms [8].

The following assumptions are introduced:

A1. The dynamic system (1) is asymptotically stable.

A2. Noises $\mathrm{<mo>\{</mo>}{\xi }_i\mathrm{<mo>\}</mo>}$ and $\mathrm{<mo>\{</mo>}{\zeta }_i\mathrm{<mo>\}</mo>}$ are statistically independent sequences with $\text{E}\mathrm{<mo>\{</mo>}{\xi }_i\mathrm{<mo>\}</mo><mo>=</mo>0,}$ $\text{E}\mathrm{<mo>\{</mo>}{\zeta }_i\mathrm{<mo>\}</mo><mo>=</mo>0,}$ $\text{E}\left\{\left.{\xi }_i^2\right\}\right.\mathrm{<mo>=</mo>}{\sigma }_{\xi }^2\mathrm{<mo><</mo>}\infty \mathrm{,}$ $\text{E}\left\{\left.{\zeta }_i^2\right\}\right.\mathrm{<mo>=</mo>}{\sigma }_{\zeta }^2\mathrm{<mo><</mo>}\infty$ a.s., where $\text{E}$ is the expectation operator.

A3. The output sequence $\left\{z_i\right\}\mathrm{,}$ is independent of noise sequence $\left\{{\xi }_i\right\}$. The noise sequences $\left\{{\xi }_i\right\}\mathrm{,}$ $\left\{{\zeta }_i\right\}$ are mutually independent.

In [18], the following objective function was proposed for estimating the parameters:

                                                            $\underset{b}{{\displaystyle \min }}\frac{{‖Y-Cb‖}_2^2}{1+\gamma +b^T{H}_{\xi }b}\mathrm{,}$                                                                  (2)

 where

$H_{\xi }^{\left(mk\right)}={\lim }_{i\to \infty }\frac{1}{N}{\sum }_{i=0}^{N-1}\left(\begin{array}{c}{\alpha }_m\\ j\end{array}\right)\left(\begin{array}{c}{\alpha }_k\\ j\end{array}\right)\frac{N-j}{N},m=\overline{\mathrm{1,}r},k=\overline{\mathrm{1,}r},$ 

$C\mathrm{<mo>=</mo>}\left({\phi }_1^T\dots {\phi }_N^T\right)\in {ℝ}^{r\times N}\mathrm{,}Y\mathrm{<mo>=</mo>}\left(y_1\dots y_N\right)\in {ℝ}^N\mathrm{,}b\mathrm{<mo>=</mo>}{\left(b^{\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}\dots b^{\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}}\right)}^T\in {ℝ}^r\mathrm{,}$ 

${\phi }_i\mathrm{<mo>=</mo>}\left({\Delta }^{{\alpha }_1}{y}_{i-1}\mathrm{,...,}{\Delta }^{{\alpha }_r}{y}_{i-1}\right)\in {ℝ}^{1\times r}\mathrm{,}\gamma \mathrm{<mo>=</mo>}{\sigma }_{\zeta }^2\mathrm{<mo>/</mo>}{\sigma }_{\xi }^2$

 Theorem 2.1. [18] Let the dynamic system described by Equation (1)with initial zero conditions and assumptions A1$–$A3 be introduced. Then, the estimate of the coefficients determined by expression (2) exists, is unique, and converges to the true value of the coefficients with probability 1, i.e.:

                                                                $\widehat{b}\mathrm{<mo\; stretchy="false">(</mo>}N\mathrm{<mo\; stretchy="false">)</mo>}\stackrel{a\mathrm{.}s\mathrm{.}}{{\displaystyle \underset{\text{N}\to \infty }{{\displaystyle \to }}}}{b}_0$                                                                      (3)

.

 Proof. The proof of the theorem is similar to the proof given in [20].

The minimum of function (2) can be found as a solution to the biased normal system of equations

                                                          $\left(C^{T}C-{\widehat{\sigma }}_{\xi }^2{H}_{\xi }\right)\widehat{b}\mathrm{<mo>=</mo>}{C}^{T}Y\mathrm{.}$                                                                (4)

 If the noise variance is unknown ${\sigma }_{\xi }^2$, then it is necessary to use the estimate of the additive noise variance ${\sigma }_{\xi }^2$. The variance estimate ${\widehat{\sigma }}_{\xi }^2$ can be found as the minimal generalized singular value

                                                             ${\widehat{\sigma }}_{\xi }\mathrm{<mo>=</mo>}{\sigma }_{\min }\left(\overline{C}\mathrm{,}{L}_{\xi }\right)\mathrm{,}$                                                                   (5)

 where ${\sigma }_{\min }\left(\overline{C}\mathrm{,}{L}_{\xi }\right)$ is the minimal generalized singular number of matrices $\overline{C}$ and $L_{\xi }$,

$\overline{C}\mathrm{<mo>=</mo>}\left(YC\right)$.

${\overline{H}}_{\xi }\mathrm{<mo>=</mo>}{\overline{L}}_{\xi }^T{\overline{L}}_{\xi }\mathrm{,}$ ${\overline{H}}_{\xi }\mathrm{<mo>=</mo>}\left(\begin{array}{cc}1+\gamma & 0\\ 0& H_{\xi }\end{array}\right)$.

In [17] a review of methods for parametr estimation integer-order autoregressions with additive noise is presented. One of the most accurate was the approach proposed in the article [21]. This article uses a generalization of this approach to the case of autoregressions with fractional order differences. The maximum value of the variance ${\sigma }_{\xi \max }^2$ is if the variance ${\sigma }_{\zeta }^2\mathrm{<mo>=</mo>0}$ is defined as

                                                           ${\sigma }_{\xi \max }^2\mathrm{<mo>=</mo>}{\sigma }_{\min }^2\left(\overline{C}\mathrm{,}{L}_{max}\right)$

where ${\sigma }_{\min }\left(\overline{C}\mathrm{,}{L}_{max}\right)$ is minimal generalised singular values of matrices $\overline{C}$ and $L_{max}\mathrm{,}$

${\overline{H}}_{max}\mathrm{<mo>=</mo>}{\overline{L}}_{max}^T{\overline{L}}_{max}\mathrm{,}$ ${\overline{H}}_{max}\mathrm{<mo>=</mo>}\left(\begin{array}{cc}1& 0\\ 0& H_{\xi }\end{array}\right)$.

The true value of the variance belongs to the interval ${\sigma }_{\xi }^2\in \left(0{\sigma }_{\xi \max }^2\right)$.

In [21], high-order Yule-Walker equations are used to determine the variance. However, this approach cannot be applied directly, since it is impossible to obtain a vector of instrumental shifts for equation (1).

Minimization (2) can be written as an eigenvector problem:

                                                            $\left({\overline{C}}^T\overline{C}-{\widehat{\sigma }}_{\xi }^2{\overline{H}}_{\xi }\right)\widehat{\overline{b}}=\mathrm{0,}$                                                                 (6)

 where $\overline{b}=\left(\begin{array}{l}-1\\ b\end{array}\right)$.

Equation (6) requires knowing not only the variance of the additive noise ${\sigma }_{\xi }^2$, but also the variance ${\sigma }_{\zeta }^2$. In order to eliminate the need to evaluate ${\sigma }_{\xi }^2$ and ${\sigma }_{\epsilon }^2$ simultaneously for fractional order autoregressions, we use generalized instrumental variables [22], the application of generalized instrumental variables for fractional order systems is considered in the article [23].

The vector of instrumental variables ${\psi }_i$ satisfies the equality

                                                          $\underset{N\to \infty }{lim}\left(C_{\psi }^T\overline{C}-{\sigma }_{\xi }^2{\overline{H}}_{\psi }\right)\overline{b}=\mathrm{0,}$                                                               (7)

 where

$C_{\psi }\mathrm{<mo>=</mo>}\left({\psi }_1^T\dots {\psi }_N^T\right)\in {ℝ}^{r\times N}\mathrm{,}$

${\psi }_i\mathrm{<mo>=</mo>}\left({\Delta }^{{\alpha }_1}{y}_{i-2}\mathrm{,...,}{\Delta }^{{\alpha }_r}{y}_{i-2}\right)\in {ℝ}^{1\times r}$, ${\overline{H}}_{\psi }\mathrm{<mo>=</mo>}\left(\begin{array}{cc}0& H_{\psi }\end{array}\right)$,

$H_{\psi }^{\left(mk\right)}={\lim }_{i\to \infty }\frac{1}{N}{\sum }_{i=0}^{N-1}\left(\begin{array}{c}{\alpha }_m\\ j-1\end{array}\right)\left(\begin{array}{c}{\alpha }_k\\ j\end{array}\right)\frac{N-j}{N},m=\overline{\mathrm{1,}r},k=\overline{\mathrm{1,}r},$

For a finite sample, equality (7) will not be strict, the problem of determining the variance estimate ${\widehat{\sigma }}_{\xi }^2$ can be described as a quadratic function minimization problem

                                                             $\underset{{\sigma }_{\xi }\in \left(\mathrm{0,}{\sigma }_{\xi \max }\right)}{{\displaystyle \min }}J\left({\sigma }_{\xi }\right)\mathrm{,}$                                                                   (8)

 where

$J\left({\sigma }_{\xi }\right)={\overline{b}}^T{\left(C_{\psi }^T\overline{C}-{\sigma }_{\xi }^2{\overline{H}}_{\psi }\right)}^T\left(C_{\psi }^T\overline{C}-{\sigma }_{\xi }^2{\overline{H}}_{\psi }\right)\overline{b}.$

Based on equations (4), (6) and (8), an iterative algorithm is proposed for estimating the parameters $\widehat{b}$ and the variance ${\widehat{\sigma }}_{\xi }^2$ .

 Step 1. Determine the maximum value of the variance ${\sigma }_{\xi \max }^2$ is defined as

                                                           ${\sigma }_{\xi \max }^2\mathrm{<mo>=</mo>}{\sigma }_{\min }^2\left(\overline{C}\mathrm{,}{L}_{max}\right)\mathrm{.}$

 Step 2. Start from a generic value ${\sigma }_{\xi }^2\in \left(0{\sigma }_{\xi \max }^2\right).$

 Step 3. Compute the parameter vector from equation (4)

                                                          $\left(C^{T}C-{\widehat{\sigma }}_{\xi }^2{H}_{\xi }\right)\widehat{b}\mathrm{<mo>=</mo>}{C}^{T}Y\mathrm{,}$

 Step 4. Compute the cost function (8)

$J\left({\sigma }_{\xi }\right)={\widehat{\overline{b}}}^T{\left(C_{\psi }^T\overline{C}-{\sigma }_{\xi }^2{\overline{H}}_{\psi }\right)}^T\left(C_{\psi }^T\overline{C}-{\sigma }_{\xi }^2{\overline{H}}_{\psi }\right)\widehat{\overline{b}}.$

 Step 5. Choose a new value ${\sigma }_{\xi }^2$ . The choice can be made using one of the methods of one-dimensional optimization.

 Step 6. Repeat steps 3$–$5 until the value associated with the minimum of is found.

2 Simulation results

 The proposed algorithm has been compared with ordinary least squares and the algorithm based on objective function (2) with a known noise variance ratio. The minimum (2) of the objective function can be found from the solution of the equation (4) or the augmented system of equations [24].

Test cases were compared by the following characteristics:the normalized root mean square error (NRMSE) of parameter estimation, defined as

$\delta b=\sqrt{{||\widehat{b}-b_0||}^2/{||b_0||}^2\cdot }100\%,$

and normalized root mean square error of modelling (NRMSEM), defined as

$\delta z=\sqrt{{||\widehat{z}-z||}^2/{||z||}^2\cdot }100\%.$

The results were based on 50 independent Monte-Carlo simulations.

 Example 1. The AR model is described by the equation

                                                      $z_i\mathrm{<mo>=</mo>0.45}{\Delta }^{-0.1}{z}_{i-1}+{\zeta }_i\mathrm{,}{y}_i\mathrm{<mo>=</mo>}{z}_i+{\xi }_i\mathrm{,}$                                                            (1)

 Noise standard deviation ratio

                                                           ${\sigma }_{\xi }\mathrm{<mo>/</mo>}{\sigma }_z\mathrm{<mo>=</mo>0.5,}\gamma \mathrm{<mo>=</mo>2.605}$

The number of data points N in each simulation was 10000.

Table 2.1 shows the mean values of tNRMSE and NRMSEM and their standard deviations.

Table 2.1

Mean values of NRMSE and NRMSEM and their standard deviations

Таблица 2.1

Средние значения NRMSE и NRMSEM и их стандартные отклонения 

Example 2. The AR model is described by the equation

                                                       $z_i\mathrm{<mo>=</mo>0.5}{\Delta }^{0.7}{z}_{i-1}+{\zeta }_i\mathrm{,}{y}_i\mathrm{<mo>=</mo>}{z}_i+{\xi }_i\mathrm{,}$                                                             (2)

 Noise standard deviation ratio

                                                           ${\sigma }_{\xi }\mathrm{<mo>/</mo>}{\sigma }_z\mathrm{<mo>=</mo>0.5,}\gamma \mathrm{<mo>=</mo>2.11}$

The number of data points N in each simulation was 2000.

Table 2.2 shows the mean values of tNRMSE and NRMSEM and their standard deviations.

Table 2.2

 Mean values of NRMSE and NRMSEM and their standard deviation

Таблица 2.2 

 Средние значения NRMSE и NRMSEM и их стандартное отклонение 

Conclusion

This paper proposed an estimation method of the parameters of fractional AR models with additive noise. The simulation results showed that the parameter estimates obtained using the proposed algorithm are highly accurate.

Further development of the proposed approach is the study of the best choice of instrumental variables and the choice of the weighting matrix.