Hybrid Model for Time Series of Complex Structure with ARIMA Components

: A hybrid model for the time series of complex structure (HMTS) was proposed. It is based on the combination of function expansions in a wavelet series with ARIMA models. HMTS has regular and anomalous components. The time series components, obtained after expansion, have a simpler structure that makes it possible to identify the ARIMA model if the components are stationary. This allows us to obtain a more accurate ARIMA model for a time series of complicated structure and to extend the area for application. To identify the HMTS anomalous component, threshold functions are applied. This paper describes a technique to identify HMTS and proposes operations to detect anomalies. With the example of an ionospheric parameter time series, we show the HMTS efﬁciency, describe the results and their application in detecting ionospheric anomalies. The HMTS was compared with the nonlinear autoregression neural network NARX, which conﬁrmed HMTS efﬁciency.


Introduction
Time series modeling and analysis are important bases for the methods of studying the processes and phenomena of different natures. They are used in various spheres of human activity (physics, biology, medicine, economics, etc.). Methods of data modeling and analysis aimed at detecting and identifying anomalies are of special actuality. The examples are the problems of the recognition of anomalies in geophysical monitoring data, such as the detection of magnetic and ionospheric storms [1][2][3][4], earthquakes [5,6], tsunamis [7,8], geological anomalies [9] and other catastrophic natural phenomena. The need to detect anomalies often arises in the medical field, for example, to detect and to identify clinical conditions of patients [10]. An important property of such methods is their ability to adapt, providing the possibility to detect and identify rapid changes in the system or object state, indicating anomaly occurrences.
As a rule, time series of empirical data have a complex non-stationary structure and contain local features of various forms. The methods for the time series analysis include deterministic [11], stochastic [12][13][14] approaches and their various combinations [15][16][17][18][19]. Traditional methods for data time series modeling and analysis (AR models, ARMA [20,21], exponential smoothing [22], stochastic approximation [13], etc.) do not allow us to describe the time series of complex structure adequately [23]. At present, hybrid approaches [16,17,19,[23][24][25][26][27][28] are widely applied. They make it possible to improve the efficiency of the procedure of data analysis in case of its complicated structure. For example, in [19], on the basis of wavelet decomposition, a technique was developed to estimate the coefficients of turbulent diffusion and power exponents from single Lagrangian trajectories of particles. Wavelet transform is a flexible tool and was applied in the paper [29] to study the relationship between vegetation and climate in India. The 2D empirical wavelet filters

Description of the Method
The time series of a complex structure may be represented as where A REG (t) = ∑ µ=1,T α µ (t) is a regular component, which is a linear combination of the components α µ (t), µ is the component number; U(t) is the anomalous component including local features of various forms occurring at random times, e(t) is the noise component, t is time.

Wavelet Series Expansion and Determination of the Model Regular Components
It is assumed that f ∈ L 2 (R)(L 2 (R) is Lebesgue space) there is a unique representation [36] f (t) = . . . + g −1 (t) + g 0 (t) + g 1 (t) + . . . , where g j ∈ W j , j ∈ Z (Z is the set of integers), g j (t) = ∑ k d j,k Ψ j,k (t), Ψ j,k = Ψ j,k k∈Z is the basis of the space W j , the coefficients d j,k = f , Ψ j,k , Ψ j,k = 2 j/2 Ψ 2 j t − k are considered as a result of mapping f into the space W j with resolution j. If Ψ ∈ L 2 (R) is R-function and the sequence Ψ j,k is a Riesz basis [37] in L 2 (R), space L 2 (R) expansion structure generated by the wavelet Ψ ∈ L 2 (R) is where W j := clos L 2 (R) Ψ j,k ; k ∈ Z , the dots above the summation sign and above the plus signs denote the direct sum.
Using expansion (2), we obtain a sequence of nested and closed subspaces V j ∈ L 2 (R), j ∈ Z defined as V j = . . .
where the space V j = clos L 2 (R) φ 2 j t − k , φ is the scaling function. Based on (2) and (3), we obtain space L 2 (R) expansion: in case of an orthogonal wavelet Ψ, we have where ⊕ is the orthogonal sum. Considering the space V j = clos L 2 (R) φ 2 j t − k with j = 0 as the base space f , and using (4) m times, we obtain the following expansion [36]: In this case, for f 0 we have the following representation: Note that, when the scaling function φ has L zero moments, i.e., 1, L and f ∈ C L (C L is the space of functions continuously differentiable by L times), then for t near 2 m k [38]: It follows from (6) that the component f −m ∈ V −m gives approximation f with resolution 2 m (it approximates the trend). The detailing component g j has the resolution Mathematics 2021, 9, 1122 4 of 18 2 −j , and approximates the local features of the scale j. Figure 1 shows the amplitudefrequency characteristics (AFC) of the scaling function (solid line) and the wavelet (dashed line) for different m, obtained for the 3rd-order Daubechies wavelet.
gives approximation f with olution 2 m (it approximates the trend). The detailing component j g has the resolu 2 j − , and approximates the local features of the scale j . Figure 1 shows the amplitu frequency characteristics (AFC) of the scaling function (solid line) and the wavelet (das line) for different m , obtained for the 3rd-order Daubechies wavelet.   [20]. Taking into account the fact that the f resolu decreases with the m increase, we define r m sequentially: Thus, we can obtain different representations of f 0 in the form (5) for different m. Obviously, it is necessary to determine the level of expansion m r , for which the component f −m r is regular. It is natural to assume that the component f −m is regular if it is strictly stationary. In this case, the problem of determining regular components is reduced to the problem of obtaining representation (5) for which the component f −m is strictly stationary. The condition of stationarity of the component f −m will allow us to identify the ARIMA model for it. Following the theory by Box and Jenkins [20], a time series is strictly stationary if its autocorrelation function (ACF) damps rapidly during average and large delays. To determine the model type (AR, MA, ARMA) and the order, ACF and partial ACF (PACF) are studied [20]. Taking into account the fact that the f resolution decreases with the m increase, we define m r sequentially: The components f −m r and g j r obtained on the basis of Algorithm 1 describe the regular changes of the time series. Then from (1) and (5), we have the representation: where , and we assume that f −m r (t) = α 1 (t), g j r (t) = α µ (t), µ = 2, T, T is the number of regular components; P j = j = −1, −(m r − 1) j = j r .

1.
We map (5) for the expansion level m = 1 for We check the condition of strict stationarity for the component f −m by estimating the numerical characteristics (analysis of ACF and PACF [20]); 3.
In the case of strict stationarity of the component f −m , we assume that it describes regular data changes (m = m r ) and go to step 5, otherwise go to step 4; 4.
If m < M, where M is the maximum level of expansion: M ≤ log 2 N (N is the time series length), we increase the expansion level by 1: m = m + 1 and return to step 2; if, m ≥ M we terminate the algorithm execution; 5.
We check the condition of strict stationarity for the detailing components g j (t) = ∑ k d j,k Ψ j,k (t), j = −1, −m r by estimating the numerical characteristics (analysis of ACF and PACF [20]). If the condition of strict stationarity is satisfied for the component g j , we take j = j r and assume that the component g j r is regular.

Estimation of the Parameters for the Model Regular Component
The components f −m r and g j r are strictly stationary, thus, we can estimate ARIMA models of order (p, ν, h) for them [20]. Then for the component for brevity, we omit index r and obtain where ω −m,k = ∇ ν c −m,k , ∇ ν is the difference operator of order ν; p, γ m 1 , . . . , γ m p are the order and the parameters of autoregression, respectively; h, θ m 1 , . . . , θ m h are the order and parameters of the moving average, respectively; a −m,k are residual errors.
In a similar way, for the component g j r (t) = ∑ k d j r ,k Ψ j r ,k (t) we omit index r and obtain where ω j,k = ∇ ν j d j,k , ∇ ν j is the difference operator of order ν j , z,, γ j 1 , . . . , γ j z are the order and the parameters of autoregression, respectively; u,, θ j 1 , . . . , θ j u are the order and the parameters of the moving average, respectively; a j,k are residual errors. From (7) to (9) we obtain the representation: where s The identification of the ARIMA model for the µ-th component requires the determination of the different order ν µ and the identification of the resulting ARMA process (model order and parameter estimation). The ARIMA model identification is described in detail in [20] and is not presented in the paper.
The diagnostic verification of each of the components f −m r and g j r models can be based on the analysis of the model residual errors. Commonly used tests based on the Mathematics 2021, 9, 1122 6 of 18 analysis of model residual errors are the cumulative fitting criterion [20] and the cumulative periodogram test [20].

Anomalous Component of the Model
The anomalous component U(t) of model (1) includes local features of various shapes occurring at random times. Therefore, the application of the parametric approach to identify it is ineffective.

Application of Threshold Functions
In the case of a nonparametric approach, following the results of [37], the function U can be approximated by threshold functions: In this case, from (7) and (10), we obtain the hybrid model of time series (HMTS) It was shown in [37] that the mappings (11) allow us to obtain approximations close to optimal ones (by minimizing the minimax risk) for a complex structure function. Moreover, the equivalence of discrete and continuous wavelet expansions [36,38] provides the opportunity to analyze a function on any resolution. In its turn, the increase in the amplitudes of the wavelet coefficients d j,k in the vicinity of local features of a function (Jaffard's theorem [39]) will provide, based on (11), their mapping into the component U of model (12).
Obviously, by applying different orthogonal wavelets Ψ we can obtain different representations (12).
We should note that due to the random nature of U, application of any thresholds T j (see (11)) is inevitably associated with erroneous decisions. In this case, the thresholds can be chosen by minimizing the posteriori risk [40].
The threshold divides the F value space of the function under analysis into two nonintersecting domains F 1 and F 2 determining anomalous and non-anomalous states, respectively. For the specific state h b , the loss average can be estimated as [40] R where ∏ bz is the loss function, P{ f ∈ F z |h b } is the conditional probability of falling within the domain F z if the state h b actually exists, b = z, b, z are the state indices ("|" denotes conditional probability).
Averaging the conditional function of the risk over all the states h b we obtain the average risk where p b is a priori probability of the state h b . If we do not know priori probabilities of the states p b , then having statistical (priori) data, we can determine posteriori probabilities P{h b | f }, b = 1, 2. Then, applying a simple loss function from (13) and (14), a posteriori risk equals

Analysis of the Model's Regular Component Errors and Detection of Anomalies
Obviously, during anomalous periods, the residual errors of the model regular component A REG (see (10)) increase. Then anomaly detection can be based on the conditional test where q ≥ 1 is the data lead step, a µ j,k are the residual errors of the µ-th component model, Q µ is the data lead length.
We can estimate the confidence interval of the predicted data [20], which is why it is logical to define the thresholds H µ as It is also possible to use the following probability limits: where u ε/2 is the quantile of the level (1 − ε/2) of standard normal distribution.

Modeling of Ionospheric Parameter Time Series
The ionosphere is the upper region of the earth's atmosphere. It is located at heights from 70 to 1000 km and higher, and affects radio wave propagation [41]. Ionospheric anomalies occur during extreme solar events (solar flares and particle ejections) and magnetic storms. They cause serious malfunctions in the operation of modern ground and space technical equipment [42]. An important parameter characterizing the state of the ionosphere is the critical frequency of the ionospheric F2-layer (foF2). The foF2 time series have a complex structure and contain seasonal and diurnal components, as well as local features of various shapes and durations occurring during ionospheric anomalies. Intense ionospheric anomalies can cause failures in the operation of technical systems. Therefore, their timely detection is an important applied problem.
In the experiments, we used hourly (1969-2019) and 15-min (2015-2019) foF2 data obtained by the method of vertical radiosonding of the ionosphere at Paratunka station (53.0 • N and 158.7 • E, Kamchatka, Russia, IKIR FEB RAS). The proposed HMTS was identified separately for foF2 hourly and 15-min data.
To identify HMTS regular components, we used the foF2 data recorded during the periods of absence of ionospheric anomalies. The application of Algorithm 1 showed that the components f −3 and g −3 are stationary (having damping ACF), thus ARIMA models can be identified for them. Figures 2 and 3 show ACF and PACF of foF2 initial time series, as well as the components f −3 and g −3 . The results confirm stationarity of the components f −3 and g −3 . An analysis of PACF shows the possibility to identify the AR models of orders 2 and 3 for the first differences of these components. The results in Figures 2 and 3 also illustrate that foF2 initial time series are non-stationary and, therefore, it is impossible to approximate them by ARIMA model without wavelet decomposition operation.   According to ratio (10) and based on the PACF of the first differences of the components 3 f − and 3 g − (Figure 3e,f), we obtain the HMTS regular component   According to ratio (10) and based on the PACF of the first differences of the components 3 f − and 3 g − (Figure 3e,f), we obtain the HMTS regular component According to ratio (10) and based on the PACF of the first differences of the components f −3 and g −3 (Figure 3e,f), we obtain the HMTS regular component Estimated parameters for s 1 −3,k and s 2 −3,k are presented in Table 1. The parameters were estimated separately for different seasons and different levels of solar activity.  Based on the data from Table 1 we obtain (1) for wintertime: for summertime and high solar activity: for summertime and low solar activity: Figure 4 shows the modeling results for HMTS regular components ( f −3 and g −3 ) during the absence of ionospheric anomalies. The model errors do not exceed the confidence interval that indicates their adequacy.   Based on the data from Table 1 we obtain (1) for wintertime: (2) for summertime and high solar activity: (3) for summertime and low solar activity:    Tables 2 and 3, and Figure 5 show the results of validation tests for the obtained models. The tests were carried out for the foF2 data that were not used at the stage of model identification. In order to verify the models, we used the cumulative fitting criterion (Tables 2 and 3), analysis of model residual error ACF (Figure 5a,b) and normalized cumulative periodogram (Figure 5c,d).        Based on the cumulative fitting criterion [20], the fitted model is satisfactory if is distributed approximately as χ 2 (Z − p − h), where Z are the considered first autocorrelations of model errors, p is the AR model order, h is the MA model order, y z (a) are the autocorrelations of model error series, n = N − ν, N is the series length, ν is the model difference order.
According to the criterion, if the model is inadequate, the average Y grows. Consequently, the model adequacy can be verified by comparing Y with the table of χ 2 distribution. The results in Tables 2 and 3 show that the Y values of the estimated models, at a significance level α = 0.05, do not exceed the table χ 2 values. The model adequacy is also confirmed by the analysis of residual error ACF (Figure 5a,b) and the normalized cumulative periodogram (Figure 5c,d). Figure 6a,b shows the results of modeling of the hourly foF2 data during the magnetic storm on 18 and 19 December 2019. Figure 6c shows the geomagnetic activity index K (K-index), which characterizes geomagnetic disturbance intensity. The K-index represents the values from 0 to 9, estimated for the three-hour interval. It is known that during increased geomagnetic activity (K > 3), anomalous changes are observed in ionospheric parameters [43]. The analysis of the results in Figure 6 shows an increase in the model errors during the increase in K-index and magnetic storm occurrence (Figure 6b). This indicates ionospheric anomaly occurrences. The results show that the HMTS allows us to detect ionospheric anomalies successfully.    Figure 7 shows the results of the application of operation (11) to 15-min foF2 data during the same magnetic storm. Based on operation (11), ionospheric anomaly occurrences are determined by the threshold function P j d j,k with the thresholds T j .
In this paper, we used the thresholds where the coefficient V = 2.3 was estimated by minimizing a posteriori risk (ratio (15)), d j,k is the average value calculated in a moving time window with the length Φ = 480 (it corresponds to the interval of 5 days).
shown in red, negative ones are shown in blue. Figure 7d shows the K-index values. The results show the occurrence of a negative ionospheric anomaly during the initial and the main phases of the magnetic storm (18 December 2019), and a positive ionospheric anomaly during the recovery phase of the storm (19 December 2019). The observed dynamics of the ionospheric parameters are characteristic of the periods of magnetic storms [43]. The results show the efficiency of HMTS application for detecting ionospheric anomalies of different intensities.

Comparison of HMTS with NARX Neural Network
To evaluate the HMTS efficiency, we compared it with the NARX neural network [44]. The NARX network is a non-linear autoregressive neural network, and it is often Positive (P j d j,k > 0) and negative (P j d j,k < 0) anomalies were considered separately. Positive anomalies (shown in red in Figure 7b) characterize the anomalous increase in foF2 values. Negative anomalies (shown in blue in Figure 7b) characterize anomalous decrease in foF2 values. To evaluate the intensity of ionospheric anomalies we used the value Assessment of the intensity of positive I + k (P j d j,k > 0) and negative I − k (P j d j,k < 0) ionospheric anomalies is shown in Figure 7c, positive anomalies are shown in red, negative ones are shown in blue. Figure 7d shows the K-index values. The results show the occurrence of a negative ionospheric anomaly during the initial and the main phases of the magnetic storm (18 December 2019), and a positive ionospheric anomaly during the recovery phase of the storm (19 December 2019). The observed dynamics of the ionospheric parameters are characteristic of the periods of magnetic storms [43]. The results show the efficiency of HMTS application for detecting ionospheric anomalies of different intensities.

Comparison of HMTS with NARX Neural Network
To evaluate the HMTS efficiency, we compared it with the NARX neural network [44]. The NARX network is a non-linear autoregressive neural network, and it is often used to forecast time series [44][45][46][47]. The architectural structure of recurrent neural networks can take different forms. There are NARX with a Series-Parallel Architecture (NARX SPA) and NARX with a Parallel Architecture (NARX PA) [44,45].
The dynamics of the NARX SPA model is described by the equation where F(·) is the neural network display function, y(k + 1) is the neural network output, x(k), x(k − 1), . . . , x(k − l x ) are neural network inputs, y(k), y(k − 1), . . . , y k − l y are past values of the time series. In NARX PA, the network input takes the network outputsŷ i =ŷ(i) instead of the past values of the time series y i = y(i), i = k, k − l y .
The neural networks were trained separately for different seasons and different levels of solar activity. During the training, we used the data for the periods without ionospheric anomalies. We obtained the networks with delays l x = l y = 2 and l x = l y = 5 for each season. The results of the networks are shown in Figure 8. Table 4 shows the standard deviations of errors (SD) of networks, which were determined as where ( ) F ⋅ is the neural network display function, ( ) The neural networks were trained separately for different seasons and different levels of solar activity. During the training, we used the data for the periods without ionospheric anomalies. We obtained the networks with delays 2 x y l l = = and 5 x y l l = = for each season. The results of the networks are shown in Figure 8. Table 4 shows the standard deviations of errors (SD) of networks, which were determined as The analysis of the results (Figure 8, Table 4) shows that the NARX SPA predicts the data with fewer errors than the NARX PA. Sending the past time series values to the NARX SPA network input (rather than network outputs) made it possible to obtain a more accurate data prediction. The comparison results of the NARX SPA with the HMTS are presented below.   The analysis of the results (Figure 8, Table 4) shows that the NARX SPA predicts the data with fewer errors than the NARX PA. Sending the past time series values to the NARX SPA network input (rather than network outputs) made it possible to obtain a more accurate data prediction. The comparison results of the NARX SPA with the HMTS are presented below. Figure 9 shows the results of ionospheric data modeling based on HMTS and NARX SPA during the periods of absence of ionospheric anomalies. The results show that the model errors have similar values for the winter and summer seasons, and vary within the interval of [−1,1], both for HMTS and NARX SPA.  Figure 9 shows the results of ionospheric data modeling based on HMTS and NARX SPA during the periods of absence of ionospheric anomalies. The results show that the model errors have similar values for the winter and summer seasons, and vary within the interval of [-1,1], both for HMTS and NARX SPA.     However, an increase in NARX SPA errors is also observed in wintertime on the eve and after the magnetic storm (Figure 10c,d). This shows the presence of false alarms. The results of detecting ionospheric anomalies based on HMTS and NARX SPA are shown in Tables 5-8. The estimates were based on statistical modeling. The HMTS results are shown for the 90% confidence interval. The analysis of the results shows that NARX The results of detecting ionospheric anomalies based on HMTS and NARX SPA are shown in Tables 5-8. The estimates were based on statistical modeling. The HMTS results are shown for the 90% confidence interval. The analysis of the results shows that NARX SPA efficiency exceeds that for HMTS during high solar activity. However, the frequency of false alarms for HMTS is significantly less than that for NARX SPA.

Conclusions
The paper proposes a hybrid model of time series of complex structure. The model is based on the combination of function expansions in a wavelet series with ARIMA models. Ionospheric critical frequency data were used to estimate the HMTS efficiency. The estimates showed:

1.
The HMTS regular component adequately describes ionospheric parameter time series during the periods without ionospheric anomalies. Application of wavelet decomposition allows us to detect regular components of ionospheric parameter time series and to use the ARIMA model; 2.
Analysis of HMTS regular component errors allows us to detect ionospheric anomalies during a magnetic storm; 3.
The HMTS anomalous component allows us to detect ionospheric anomalies of different intensities by threshold functions.
Comparison of HMTS with NARX with Series-Parallel Architecture confirmed the HMTS efficiency to detect anomalies in the ionospheric critical frequency data. The results of the experiments showed that the efficiency of the NARX neural network slightly exceeds that of HMTS (about 2-3%) during high solar activity. However, the frequency of false alarms in NARX is significantly higher (about 15%). During the periods of low solar activity, the efficiency of HMTS exceeds that of NARX.
The HMTS can be used for modeling and analysis of time series of complex structure, including seasonal components and local features of various forms.