Accelerated Maximum Entropy Method for Time Series Models Estimation

Abstract

The work is devoted to the development of a maximum entropy estimation method with soft randomization for restoring the parameters of probabilistic mathematical models from the available observations. Soft randomization refers to the technique of adding regularization to the functional of information entropy in order to simplify the optimization problem and speed up the learning process compared to the classical maximum entropy method. Entropic estimation makes it possible to restore probability distribution functions for model parameters without introducing additional assumptions about the likelihood function; thus, this estimation method can be used in problems with an unspecified type of measurement noise, such as analysis and forecasting of time series.

1. Introduction

This work is devoted to the development of a maximum entropy estimation method with soft randomization for restoring the values of the parameters of probabilistic mathematical models from the available observations. Soft randomization is understood as the technique of adding regularization to the information entropy functional in order to simplify the optimization problem and speed up the learning process was originally proposed in a recent work [1].

When solving regression problems, it is generally accepted to separate the classical statistical approach, in which the model parameters are considered deterministic, and their values are estimated, for example, using the least squares method (LS) or the maximum likelihood method (ML) [2] and the probabilistic approach, according to which the model parameters are assumed to be random variables, and the result of estimation is their probabilistic characteristics [3]. This approach corresponds, for example, to Bayesian estimation and entropy-robust estimation—the general maximum entropy method (GME) [4]. One of the differences between the statistical and probabilistic approaches is that the latter allows the use of probabilistic characteristics and properties of random variables to refine estimates and subsequent forecasts using them.

In practice, the transition to probabilistic estimation is most common in problems where the prerequisites of the classical least squares method and the Gauss–Markov theorem [5] are violated. Examples are problems with a small amount of input data, when it becomes difficult to test statistical hypotheses, and problems with non-standard types of errors or, in general, unstructured noise. Estimates obtained by probabilistic methods for such problems may turn out to be biased, but more efficient than classical methods estimates, such as LS, in terms of the mean square error of estimates [6].

It is worth noting that entropy and other entropy-based metrics are widely used in the field of data analysis, in general, and time series analysis, in particular. For instance, in works [7,8,9,10], the entropy concepts are used as measures of time series dependence as an alternative to correlation-based metrics.

In [11,12], the concept of soft entropy estimation was considered, including obtaining entropy-optimal probability density functions for model parameters in general form, and the effectiveness of the proposed method of approximate maximum entropy estimation was also investigated. This method was compared with the maximum likelihood method and Bayesian estimation method, both experimentally and theoretically for special cases. Estimation methods were tested on the example of a linear regression problem with various types of errors. The experiments carried out make it possible to distinguish a class of problems in which the proposed method is more efficient than its analogs. These problems include cases where the distribution of model errors differs from the standard normal distribution in terms of skewness and kurtosis. In addition, the efficiency of maximum entropy estimation becomes more noticeable under the conditions of a small amount of input data.

This paper presents the results of a study of the practical applicability of the soft maximum entropy method in the problems of analysis and forecasting of time series. The aim of this study is to analyze the effectiveness of the proposed approach for solving practical problems of time series analysis.

The rest of the article is organized as follows: in Section 2 and Section 3, we provide a brief time series analysis review and revisit ARMA model estimation methods; in Section 4, we present a variation of traditional maximum entropy method and apply it to the ARMA model estimation problem; Section 5 and Section 6 contain the results of experiments with synthetic data and real data, respectively; and Section 7 is the conclusion and final summary of the presented work.

2. Time Series Analysis Review

One of the most versatile and widely used models for time series analysis is the autoregressive moving average, or ARMA model for short [13,14]. This model is used to build a functional dependence of some target variable based on its past values and on the values of the noise of the model.

There is also a modification of the classical ARMA model with the addition of an integrated component of the series (Autoregressive Integrated Moving Average—ARIMA) in order to model non-stationary series by calculating the differences in a certain order from its original values [15,16]. To restore dependence on other variables, terms for the regressors of the model (ARIMAX) are added to the model, traditionally denoted by the X matrix. In other words, the construction of the most universal time series models is based on the classical ARMA model.

A second widely known group of models for time series analysis is the autoregressive conditional heteroscedasticity (ARCH) models and generalized ARCH (GARCH) models [17,18]. This type of model is used to restore the dependence of the conditional variance of a series on its past values. That is, they are rather models for analyzing the type and distribution of the errors in a series, and, therefore, are often used in conjunction with models of the series itself, such as ARMA/ARIMA [19].

In this paper, we propose a new estimation method for restoring the model parameters of some functional dependence of the time series; therefore, in the future, we will consider the ARMA model as one of the most universal for time series analysis.

3. Building an ARMA Model and Estimation Methods

The ARMA model is a combination of simpler models, namely, the autoregressive model (AR) and the moving average (MA) model, which, in general, have the following form:

$$\begin{array}{cc}\hfill AR\left(p\right):y_t=& a_0+a_1{y}_{t-1}+a_2{y}_{t-2}+\cdots +a_p{y}_{t-p}+{\epsilon }_t\hfill \\ \hfill MA\left(q\right):y_t=& b_0+b_1{\epsilon }_{t-1}+b_2{\epsilon }_{t-2}+\cdots +b_q{\epsilon }_{t-q}+{\epsilon }_t.\hfill \end{array}$$

Correspondingly:

$$ARMA(p,q):y_t=a_0+a_1{y}_{t-1}+\cdots +a_p{y}_{t-p}+b_1{\epsilon }_{t-1}+b_2{\epsilon }_{t-2}+\cdots +b_q{\epsilon }_{t-q}+{\epsilon }_t.$$

Here, p and q are the hyperparameters of the model. The orders of the model are chosen, as a rule, by analyzing the plots of the autocorrelation function (ACF) and partial autocorrelation function (PCF) [20]. Another way to choose the order of the model is to calculate the information criteria Akaike’s information criterion (AIC) and Bayesian information criterion (BIC) [21]. The preferred models are those that exhibit lower values of these criteria.

We restrict ourselves here to this brief description of the methods for choosing the model order without detailed formulas for autocorrelation functions and information criteria, which is due to the fact that the object of study in this paper is the method of estimating the model parameter, and not its order or structure.

The estimation of the model parameters is understood as the restoration of the values of the parameters $a_0,a_1,\cdots ,a_p$ and $b_1,\cdots ,b_q$ from a sample of available observations. Unlike classical linear regression, the traditional least squares method becomes inapplicable for the ARMA model due to the recursive structure of the model itself.

The most common method for estimating the parameters of an ARMA model is the maximum likelihood method (maximum likelihood estimation—MLE); this method is the “standard” for building ARMA models in most development environments, such as R, Python, MATLAB, etc. [22,23]. The method uses the likelihood or log-likelihood function and computational algorithms to find the maximum point. Despite the widespread use of the MLE method, it has two significant drawbacks that make it difficult to apply in practice.

The first disadvantage is the need to make assumptions about the true form of the likelihood function for the model noise. The traditional choice here is the normal distribution; however, there are problems in time series analysis in which the errors are distributed differently, or the size of the input data sample is insufficient to build and test such hypotheses [24,25]. This feature is characteristic of all estimation methods that use a predetermined form of the likelihood function for optimization, that is, for Bayesian estimation too [26].

Another disadvantage of the MLE method is the point nature of the resulting estimates. As a result of solving the optimization problem, one single point will be selected; this point will set the vector of model parameters and determine one single curve along which the process can develop. On the other hand, in Bayesian estimation, the result of the estimate is the distribution functions for the model parameters, which allows building confidence intervals in order to refine the forecasts.

Among the alternative methods for estimating the parameters of ARMA models, one can distinguish the method of moments based on the use of the values of the ACF/PCF characteristics and sample moments [27], the method based on the Yule–Walker estimation equations, and the analytical derivation of a solution for the model parameters [28], as well as the so-called innovations algorithm based on successive recursive approximation of the solution [29]. Some of the most recent works in the field of robust parameter estimation for time series models are [30,31].

Another alternative for searching for the parameters of the ARMA model can be the method of soft maximum entropy estimation proposed in this study. This method is an approximation of the classical maximum entropy method [32,33] in order to simplify the optimization problem and speed up calculations [1].

4. Method of Soft Maximum Entropy Estimation for ARMA Model

The method of approximate maximum entropy estimation proposed in this paper is formulated for the regression recovery problem of the following type:

$$y=F\left(X,a\right)+\epsilon ,a\in A,$$

where F is a vector function and $\epsilon \in E$ are measurement noise.

The traditional maximum entropy estimation method is formulated as follows [4]:

$$H\left[P\left(a\right),Q\left(\epsilon \right)\right]=-{\int }_{A}P\left(a\right)lnP\left(a\right)da-{\int }_{E}Q\left(\epsilon \right)lnQ\left(\epsilon \right)d\epsilon \Rightarrow \underset{P,Q}{max},$$

$${\int }_{A}P\left(a\right)da=1,{\int }_{E}Q\left(\epsilon \right)d\epsilon =1,$$

under the constraints to the average model output (first-moment balance):

$${\mathbb{E}}^P\left[F_i\left(X,a\right)\right]+{\mathbb{E}}^Q\left[{\epsilon }_i\right]=y_i,i=1,\cdots ,n.$$

The optimization problem listed above has the solution in the general form using Lagrange multipliers the search for which is the main computational task carried out by the substitution of the solution for P and Q into the balance constraints, i.e., the system of integral equations.

The accelerated variation of the maximum entropy estimation method [1] consists of optimizing the modified information entropy functional with the addition of regularization terms:

$$J\left[P\left(a\right),Q\left(\epsilon \right)\right]=H\left[P\left(a\right),Q\left(\epsilon \right)\right]-{\overline{N}}_a\left[P\left(a\right)\right]-{\overline{N}}_{\epsilon }\left[Q\left(\epsilon \right)\right]\Rightarrow \underset{P,Q}{max}$$

$${\int }_{A}P\left(a\right)da=1,{\int }_{E}Q\left(\epsilon \right)d\epsilon =1.$$

Here, ${\overline{N}}_a\left[P\left(a\right)\right]$ and ${\overline{N}}_{\epsilon }\left[Q\left(\epsilon \right)\right]$ are the average values for the Hölder model output norm and noise norm:

$$N_a\left(a\right)={∥F\left(X,a\right)-y∥}_H,N_{\epsilon }\left(\epsilon \right)={∥\epsilon ∥}_H.$$

That is:

$$\begin{array}{cc}\hfill {\overline{N}}_a\left[P\left(a\right)\right]& =\mathbb{E}\left[N_a\left(a\right)\right]={\int }_{A}P\left(a\right)N_a\left(a\right)da,\hfill \\ \hfill {\overline{N}}_{\epsilon }\left[Q\left(\epsilon \right)\right]& =\mathbb{E}\left[N_{\epsilon }\left(\epsilon \right)\right]={\int }_{E}Q\left(\epsilon \right)N_{\epsilon }\left(\epsilon \right)d\epsilon .\hfill \end{array}$$

Thus, the idea of so-called soft maximum entropy estimation is to calculate the approximate maximum of the entropy functional with the minimization of the average model error norm instead of solving the system of balance constraints.

In [11,12], the solution of the optimization problem for the soft maximum entropy estimation method was obtained in the general form:

$$\begin{array}{cc}\hfill P^{*}\left(a\right)& =\frac{exp\left(-N_a\left(a\right)\right)}{{\int }_{A}exp\left(-N_a\left(a\right)\right)da},\hfill \\ \hfill Q^{*}\left(\epsilon \right)& =\frac{exp\left(-N_{\epsilon }\left(\epsilon \right)\right)}{{\int }_{E}exp\left(-N_{\epsilon }\left(\epsilon \right)\right)d\epsilon }.\hfill \end{array}$$

The function $P^{*}\left(a\right)$ can be treated as a posterior distribution function of the parameters and, similar to the Bayesian approach, can be used to calculate the average values of the parameters or MAP estimates, while additional assumptions regarding the error likelihood function were not used to obtain them.

A feature of the proposed estimation method is the presence of a direct functional dependence, $y=F\left(X,a\right)+\epsilon$. If the autoregressive component of the time series model corresponds to this structure and is observable, then the presence of the lag terms of the moving average greatly complicates such a representation, as a result of which the general maximum entropy estimation method becomes inapplicable for the ARMA model explicitly.

Let us consider the representation of the ARMA model through the lag operator, L:

$$L{x}_t=x_{t-1},L^2{x}_t=x_{t-2},\cdots ,L^k{x}_t=x_{t-k}.$$

Then, omitting (without loss of generality) the constant, the autoregressive models for the variable x will have the following form:

$$\begin{array}{cc}\hfill AR\left(1\right):& \left(1-aL\right)x_t={\epsilon }_t\hfill \\ \hfill AR\left(p\right):& \left(1-a_{1}L-a_2{L}^2-\cdots -a_p{L}^p\right)x_t={\epsilon }_t\hfill \\ & 1-a_{1}L-a_2{L}^2-\cdots -a_p{L}^p=a^p\left(L\right).\hfill \end{array}$$

Similarly, we transform the moving average models:

$$\begin{array}{cc}\hfill MA\left(1\right):& x_t=\left(1+bL\right){\epsilon }_t\hfill \\ \hfill MA\left(q\right):& x_t=\left(1+b_{1}L+b_2{L}^2+\cdots +b_q{L}^q\right){\epsilon }_t\hfill \\ & 1+b_{1}L+b_2{L}^2+\cdots +b_q{L}^q=b^q\left(L\right).\hfill \end{array}$$

By combining the representations for AR(p) and MA(q), we obtain:

$$ARMA(p,q):a^p\left(L\right)x_t=b^q\left(L\right){\epsilon }_t.$$

Here, we consider the important property of the invertibility of the ARMA model using the example of individual first-order models.

For example, for AR(1), the following statement is true:

$$\begin{array}{cc}\hfill AR\left(1\right):& \left(1-aL\right)x_t={\epsilon }_t\hfill \\ & x_t={\left(1-aL\right)}^{-1}{\epsilon }_t\to MA(\infty )=b\left(L\right){\epsilon }_t,\hfill \end{array}$$

where $b\left(L\right)=1+b_{1}L+b_2{L}^2+\cdots +b_k{L}^k+\cdots$, and $b_k=a^k\forall k=1,\cdots$.

The proof of this fact follows from the following expression, if we substitute the required form of the coefficients, $b_k$:

$$\begin{array}{cc}\hfill \left(1-aL\right)b\left(L\right)& =\left(1-aL\right)\left(1+b_{1}L+b_2{L}^2+\cdots +b_k{L}^k+\cdots \right)\hfill \\ & =\left(1-aL\right)\left(1+aL+a^2{L}^2+\cdots +a^k{L}^k+\cdots \right)\hfill \\ & =\left(1-aL+aL-a^2{L}^2+a^2{L}^2-\cdots \right)=1-\underset{k\to \infty }{lim}{a}^k{L}^k=1,\mathrm{while}\parallel a\parallel <1.\hfill \end{array}$$

Furthermore, similarly, for MA(1), it is true:

$$\begin{array}{cc}\hfill MA\left(1\right):& x_t=\left(1+bL\right){\epsilon }_t\hfill \\ & {\left(1+bL\right)}^{-1}{x}_t={\epsilon }_t\to AR(\infty )=a\left(L\right)x_t,\hfill \end{array}$$

where $a\left(L\right)=1-a_{1}L-a_2{L}^2-\cdots -a_k{L}^k-\cdots$, and $a_k=-{(-b)}^k\forall k=1,\cdots$.

The proof is carried out in a similar way, substituting the required form of the coefficients, $a_k$:

$$\begin{array}{cc}\hfill \left(1+bL\right)a\left(L\right)& =\left(1+bL\right)\left(1-a_{1}L-a_2{L}^2-\cdots -a_k{L}^k-\cdots \right)\hfill \\ & =\left(1+bL\right)\left(1+(-b)L+{(-b)}^2{L}^2+\cdots +{(-b)}^k{L}^k+\cdots \right)\hfill \\ & =\left(1+bL\right)\left(1-bL+b^2{L}^2-b^3{L}^3+\cdots \right)\hfill \\ & =\left(1+bL-bL+b^2{L}^2-b^2{L}^2+\cdots \right)=1+\underset{k\to \infty }{lim}{b}^k{L}^k=1,\mathrm{while}\parallel b\parallel <1.\hfill \end{array}$$

The considered statements demonstrate the invertibility property for autoregressive models when the autoregressive component can be represented as an infinite series of moving average terms and vice versa.

It turns out that the more general models AR(p) and MA(q) are also invertible under certain restrictions on the values of the coefficients [34]. In general, the condition for the invertibility of the ARMA(p,q) model is the absence of identical roots of the characteristic polynomials, $a^p\left(L\right)$ and $b^q\left(L\right)$, and roots exceeding modulo 1, similar to the considered cases for a single coefficient in AR(1) and MA(1).

Then, we can represent the general ARMA model in the following form:

$$\begin{array}{cc}\hfill ARMA(p,q):& a^p\left(L\right)x_t=b^q\left(L\right){\epsilon }_t\hfill \\ & x_t={\left(a^p\right)}^{-1}\left(L\right)b^q\left(L\right){\epsilon }_t=c\left(L\right){\epsilon }_t\hfill \\ & b^q\left(L\right)=a^p\left(L\right)c\left(L\right),\hfill \end{array}$$

or vice versa:

$$\begin{array}{cc}& {\left(b^q\right)}^{-1}\left(L\right)a^p\left(L\right)x_t={\epsilon }_t=d\left(L\right)x_t,\hfill \\ & a^p\left(L\right)=b^q\left(L\right)d\left(L\right)\hfill \end{array}$$

Thus, it is possible to represent the ARMA model in such a way that all its terms are observable, except for one single term for the residuals of the model, similar to the classical regression problem.

It is on the possibility of such a representation that the Hannan–Rissanen algorithm is based on iterative calculation of parameter estimates for the ARMA model [35]. By analogy with this algorithm, we represent the problem of maximum entropy estimation for the ARMA model as a sequence of several stages:

1.

Build an AR model of a higher order (for example, for the original ARMA(p,q) model at this stage, starting with an estimate of AR(p+q));

2.

Calculate, using the model from step 1, the values of the residuals, ${\epsilon }_t$, for each time, t;

3.

Using a sample of values of the target variable X and the found values of the residuals, build an ARMA(p,q) model:

$$x_t\sim x_{t-1},\cdots ,x_{t-p},{\epsilon }_{t-1},\cdots ,{\epsilon }_{t-q}.$$

After step 3, new estimates for the residuals of the model can be used to further improve the estimates of the parameters.

Thus, at each stage, the parameters of the model are estimated, in which all terms are observable values: first, this is the target variable itself and its past values; then, the calculated values of the residuals of the model. Accordingly, this algorithm can use the maximum entropy estimation method, with both based on the original maximum entropy method (general maximum entropy—GME) and its modification proposed in this paper—the soft maximum entropy estimation method.

5. Results of Experiments with Model Data

This section presents the results of testing the performance and efficiency of the proposed algorithm in model experiments on synthetic datasets.

Testing is carried out on several time series models: AR(1), MA(1), ARMA(1,1), and, in the general case, implementation for ARMA(p,q), which is possible. In the first series of experiments, time series data with a given model are generated under standard noise corresponding to a normal distribution, ${\epsilon }_t\sim N(\mu =0,{\sigma }^2=0.1)$, as shown in Figure 1. In the second series of experiments, the distribution of measurement errors is replaced by a non-standard and asymmetric chi-square distribution—${\epsilon }_t\sim {\chi }^2(\nu =0.1)$, with examples of noisy data that include outlier processing shown in Figure 2. Finally, the third series of experiments demonstrates the example of other non-standard errors, in this case, exponential distribution ${\epsilon }_t\sim exp(\mu =0.2)$. The purpose of the experiments is to compare the effectiveness of various estimation methods under the conditions of different non-standard measurement noises.

As can be seen in Figure 1 and Figure 2, visually noisy data are practically indistinguishable from a series with standard errors; however, this can significantly affect the results of estimating model parameters, for example, due to the use of an incorrect likelihood function in the MLE method, or under assumptions about the normality of the residuals in the OLS method.

Software development and experimental testing is carried out on a standard PC with an Intel(R) Core(TM) i7-9700F CPU 3.00 Hz and 32 GB of RAM in a MATLAB software environment, version R2019b, and Econometrics Toolbox is used to generate and evaluate ARMA models. In addition, by setting the initial position of the random number generator (rng(2022)), all the results obtained are reproducible within the standard deviations.

The results of the first series of experiments are shown in

Table 1. Estimation results for the AR(1) model.

$AR\left(1\right):y_t=\mu +a{y}_{t-1}+{\epsilon }_t,{\epsilon }_t\sim N(0,{\sigma }^2=0.1)$
	$\mu =0.3$		$a=0.5$
	Mean	MSE	Mean	MSE
MLE	0.31083	0.00211	0.48587	0.00433
SME $L_2$	0.31270	0.00218	0.48286	0.00442
SME $L_1+L_2$	0.31195	0.00252	0.48343	0.00512

,

Table 2. Estimation results for the MA(1) model.

$MA\left(1\right):y_t=\mu +b{\epsilon }_{t-1}+{\epsilon }_t,{\epsilon }_t\sim N(0,{\sigma }^2=0.1)$
	$\mu =0.3$		$b=-0.5$
	Mean	MSE	Mean	MSE
MLE	0.30120	0.00013	−0.51495	0.00442
SME $L_2$	0.30123	0.00014	−0.46826	0.00588
SME $L_1+L_2$	0.30125	0.00020	−0.46476	0.00669

and

Table 3. Estimation results for the ARMA(1,1) model.

$ARMA(1,1):y_t=\mu +a{y}_{t-1}+b{\epsilon }_{t-1}+{\epsilon }_t,{\epsilon }_t\sim N(0,{\sigma }^2=0.1)$
	$\mu =0.3$		$a=0.8$		$b=-0.2$
	Mean	MSE	Mean	MSE	Mean	MSE
MLE	0.33200	0.01036	0.78004	0.00434	−0.19619	0.00987
SME $L_2$	0.33500	0.01102	0.77804	0.00461	−0.19496	0.00993
SME $L_1+L_2$	0.33423	0.01398	0.77857	0.00585	−0.19527	0.01188

for the AR(1), MA(1), and ARMA(1,1) models, respectively. The tables for each estimation method show the mean values (mean) of the estimated parameters and their standard deviation (MSE) from the true values based on the result of a series of 500 experiments with a sample size of 200 points. The estimation methods used are designated as MLE (maximum likelihood method), SME L2 (soft maximum entropy estimation method with quadratic regularization norm), and SME L1 + L2 using a linear combination of the L1 and L2 norms.

Before proceeding to the analysis of the results given in Table 1, Table 2 and Table 3, it is worth noting that the application of the soft maximum entropy estimation method with the L2 norm is equivalent to the application of the least squares method. In addition, with normal errors, the maximum likelihood estimates are also the same as the least squares estimates. Therefore, in this series of experiments, where normal errors are used, the first two lines in Table 1, Table 2 and Table 3 practically coincide. The results obtained are consistent with the results of a similar experiment in [36], which also used 200 points of the time series and 500 repetitions for the same models with the same parameter values: the MSE value for the model parameters is about 0.004, and for the constant, it about 0.002.

As can be seen in the data in Table 1, Table 2 and Table 3, the standard MLE method turned out to be the most effective estimation method in all three models because this method shows the lowest MSE error. The method of soft maximum entropy estimation with a quadratic norm coincides with this result, and with a combined norm it is slightly inferior in efficiency in terms of MSE. This result is quite natural and expected since it is known that in the case of normal errors, it is the classical estimation methods (OLS and MLE) that are unbiased and most effective for estimating regression parameters. However, the method proposed in this paper based on approximate maximum entropy estimation shows results that are close in efficiency, which proves its efficiency for time series models.

Table 4. Estimation results for the AR(1) model under ${\chi }^2$-noise.

$AR\left(1\right):y_t=\mu +a{y}_{t-1}+{\epsilon }_t,{\epsilon }_t\sim {\chi }^2\left(0.1\right)$
	$\mu =0.3$		$a=0.5$
	Mean	MSE	Mean	MSE
MLE	0.41299	0.02388	0.31391	0.06467
SME $L_2$	0.41643	0.02442	0.30817	0.06614
SME $L_1+L_2$	0.38428	0.01564	0.35991	0.04327

,

Table 5. Estimation results for the MA(1) model under ${\chi }^2$-noise.

$MA\left(1\right):y_t=\mu +b{\epsilon }_{t-1}+{\epsilon }_t,{\epsilon }_t\sim {\chi }^2\left(0.1\right)$
	$\mu =0.3$		$b=-0.5$
	Mean	MSE	Mean	MSE
MLE	0.35323	0.00318	−0.29074	0.05228
SME $L_2$	0.35327	0.00319	−0.27237	0.05906
SME $L_1+L_2$	0.34370	0.00215	−0.27470	0.05807

and

Table 6. Estimation results for the ARMA(1,1) model under ${\chi }^2$-noise.

$ARMA(1,1):y_t=\mu +a{y}_{t-1}+b{\epsilon }_{t-1}+{\epsilon }_t,{\epsilon }_t\sim {\chi }^2\left(0.05\right)$
	$\mu =0.3$		$a=0.8$		$b=-0.2$
	Mean	MSE	Mean	MSE	Mean	MSE
MLE	0.34022	0.01316	0.77824	0.00528	−0.28223	0.02000
SME $L_2$	0.34038	0.01439	0.77811	0.00582	−0.28425	0.02164
SME $L_1+L_2$	0.33783	0.01518	0.77823	0.00617	−0.26824	0.01918

show the experimental results for the same AR(1), MA(1), and ARMA(1,1) models, but under the conditions of non-standard measurement noise corresponding to a chi-square distribution. Likewise, the results of experiments with exponential measurement noise are demonstrated

Table 7. Estimation results for the AR(1) model under exp-noise.

$AR\left(1\right):y_t=\mu +a{y}_{t-1}+{\epsilon }_t,{\epsilon }_t\sim exp\left(0.2\right)$
	$\mu =0.3$		$a=0.5$
	Mean	MSE	Mean	MSE
MLE	0.80318	0.25697	−0.00930	0.26453
SME $L_2$	0.80318	0.25691	−0.00927	0.26439
SME $L_1+L_2$	0.79710	0.25090	−0.00731	0.26251

,

Table 8. Estimation results for the MA(1) model under exp-noise.

$MA\left(1\right):y_t=\mu +b{\epsilon }_{t-1}+{\epsilon }_t,{\epsilon }_t\sim exp\left(0.2\right)$
	$\mu =0.3$		$b=-0.5$
	Mean	MSE	Mean	MSE
MLE	0.79594	0.24670	−0.00886	0.24663
SME $L_2$	0.79587	0.24665	−0.00664	0.24848
SME $L_1+L_2$	0.79129	0.24220	−0.00428	0.25104

and

Table 9. Estimation results for the ARMA(1,1) model under exp-noise.

$ARMA(1,1):y_t=\mu +a{y}_{t-1}+b{\epsilon }_{t-1}+{\epsilon }_t,{\epsilon }_t\sim exp\left(0.2\right)$
	$\mu =0.3$		$a=0.8$		$b=-0.2$
	Mean	MSE	Mean	MSE	Mean	MSE
MLE	0.38372	0.02422	0.77478	0.00643	−0.33653	0.03244
SME $L_2$	0.38376	0.02558	0.77488	0.00682	−0.33795	0.03415
SME $L_1+L_2$	0.37813	0.02683	0.77610	0.00754	−0.33177	0.03327

. In this series of experiments, the classical OLS and MLE estimation methods turn out to be less effective due to the violation of the hypothesis about the normality of regression residuals.

From the results of Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9, it can be noted that in this experiment all estimates are biased due to the asymmetry of errors; however, it is the estimates using the soft maximum entropy estimation method that show the smallest deviation from the true values, although their bias does not differ in order of magnitude from the rest of the methods. Under conditions of ${\chi }^2$-noise (Table 4, Table 5 and Table 6), all estimation methods showed a comparable decrease in efficiency compared to similar experiments in Table 1, Table 2 and Table 3, so the value of the dispersion of estimates increased by about 1.5–2 times for the autoregressive model AR(1) and for the moving average model MA(1). To a lesser extent, the combined ARMA(1,1) model is prone to errors, which is probably due to the nature of the Hannan–Rissanen estimates, where the moving average component is restored through the estimates of autoregressive models of higher orders. The estimation results in Table 7, Table 8 and Table 9 are highly biased relative to the true values due to the asymmetry and tailness of the exponential distribution with significant mean value. In this case, all the methods under consideration suffer a considerable drop in efficiency as well as estimation bias.

Next, using the asymptotic normality of the estimates by the method of soft maximum entropy, we analyze the statistical significance of the results obtained for ${\chi }^2$-noise by testing the hypothesis that the means of two normal distributions are equal. To do this, for each estimated parameter, the null hypothesis, $H_0:{\mu }_{sme}={\mu }_{MLE}$, is tested against the alternative hypothesis about the inequality of means, $H_1:{\mu }_{sme}<{\mu }_{MLE}$, using Student’s t-test. In this case, the symbols $\mu$ denote not the average values of the parameters, but their average deviations (the so-called bias) from the true values given in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6 for all tested models. The results of hypothesis testing are presented in

Table 10. Statistical significance of the estimation results for ${\chi }^2$-noise.

	SME $L_2$ vs. MLE			SME $L_1+L_2$ vs. MLE
	$\mu$	$a$	$b$	$\mu$	$a$	$b$
$AR\left(1\right)$	∘	•	–	•	•	–
$MA\left(1\right)$	•	–	∘	•	–	∘
$ARMA(1,1)$	∘	∘	∘	•	∘	•

. The solid black circle indicates cases in which the null hypothesis is rejected at a significance level of 95%; therefore, the difference in the results of the estimates is significant.

During the calculations, the statistical significance of the difference in the estimation results was confirmed not in all of the experiments and not for all of the parameters. So, for example, estimates of all three parameters of the ARMA model turned out to be indistinguishable for the likelihood method and entropy estimation with the L2 regularization norm. On the other hand, entropy estimation with the combined norm shows more stable results, and the improvement in estimates was recorded in five cases out of seven. Therefore, this method will also be applied to experiments with real data.

6. Results of Experiments with Real Data

As an example of real data, we will use the NASDAQ stock exchange composite index data collection (Data_EquityIdx dataset in MATLAB since version R2012). The time series is the consecutive daily values of the index at the close of the exchange, and there are no gaps in the data. A graph of the NASDAQ index values from January 1990 to November 1997 is shown in Figure 3a, and A corresponding plot of the first difference in values is shown in Figure 3b.

The selected segment is 2000 time reports, and we will use it to test the proposed estimation method. To do this, the sample is divided in a ratio of 80% and 20% into the training set, by which the model will be evaluated, and into the test set, by which we will calculate the average forecast error.

The model is chosen according to the Box–Jenkins methodology [37] based on the values of autocorrelation (ACF) and partial autocorrelation (PCF) functions. Graphs of the corresponding functions for the difference in index values are shown in Figure 4. Analysis of these graphs shows a probable ARMA(1,1) model or, taking into account the previously calculated difference from the original values, we obtain the ARIMA(1,1,1) model.

The model parameters are estimated in three ways: the first and second are the soft maximum entropy estimation (SME) method proposed in this paper with regularization norm $L_2$ and combined $L_1+L_2$ norm, and the third is the maximum likelihood method (MLE), which is standard for constructing ARMA models in MATLAB (arima function from Econometric Toolbox). The calculation results are shown in

Table 11. Parameter estimation results for the ARMA(1,1) model.

$ARMA(1,1):y_t=\mu +a{y}_{t-1}+b{\epsilon }_{t-1}+{\epsilon }_t$
	$\mu$	a	b	${\sigma }^2$
MLE	0.4554	−0.000878	0.1839	33.313
SME $L_2$	0.4731	−0.0308	0.2137	33.454
SME $L_1+L_2$	0.4900	−0.0230	0.2079	33.454

.

On the estimation interval, all three methods showed similar results, the values of the parameters are in the range of 5–10%, and the variance in the residuals differs slightly. Moreover, the parameter a for the first-order autoregression in all three cases turned out to be insignificant (using Student’s t-test).

The resulting estimates were used to build a forecast and compare with real values from the remaining 20% of the data sample. In Figure 5, the forecast results are presented using estimates obtained by the soft maximum entropy estimation method, and the average curve (Avg. forecast) and the 95% confidence interval tube (95% CI) are marked on the graph. The maximum point of the posterior distribution density (MAP) obtained by sampling the entropy-optimal density for the model parameters was used as a point estimate of the model parameters, and the 95 percent interval was calculated over the entire ensemble of sampled values.

As in the estimation interval, during testing, all the obtained models also showed close results; on the graph, they would be indistinguishable at a given scale. Thus, the value of the root-mean-square error RMSE for all points of the forecasting stage was 151.73 for estimates using the MLE method, and for estimates using the SME method, the values were 151.06 and 147.61 using the L2 and L1 + L2 norms, respectively.

Taking into account the true values of the NASDAQ index in this interval, we can say that the three forecasts received are within 2–3% of each other. The normalized root-mean-square error (NRMSE) values obtained by dividing by the mean of the predictor are 0.1124, 0.1119, and 0.1093 for the MLE, SME L2, and SME L1 + L2 estimates, respectively.

7. Conclusions

In this work, the concept of applying the method of soft maximum entropy estimation to build time series models, such as an autoregressive model, a moving average model, and a combined ARMA model, was developed. The main advantage of the proposed method compared to the traditional maximum entropy estimation method is the elimination of solving the system of balance constraints. The disadvantage of the proposed method is the necessity for explicit regression form with observed variables.

Experimental testing of the proposed approach was carried out, both for synthetic data and for real data on the NASDAQ exchange index. The simulation of synthetic data was carried out in the presence of measurement noise, both according to the law of the standard normal distribution, and noise of a non-standard asymmetric type, according to the chi-square law. The results of the experiments demonstrate the operability of the proposed approach and an efficiency close to the classical methods of estimation.