Robust filtering of sequences with periodically stationary multiplicative seasonal increments

We study stochastic sequences $\xi(k)$ with periodically stationary generalized multiple increments of fractional order which combines cyclostationary, multi-seasonal, integrated and fractionally integrated patterns. We solve the filtering problem for linear functionals constructed from unobserved values of a stochastic sequence $\xi(k)$ based on observations with the periodically stationary noise sequence. For sequences with known matrices of spectral densities, we obtain formulas for calculating values of the mean square errors and the spectral characteristics of the optimal estimates of the functionals. Formulas that determine the least favorable spectral densities and minimax (robust) spectral characteristics of the optimal linear estimates of the functionals are proposed in the case where spectral densities of sequences are not exactly known while some sets of admissible spectral densities are given.


Introduction
Non-stationary and long memory time series models are wildly used in different fields of economics, finance, climatology, air pollution, signal processing etc. (see, for example, papers by Dudek and Hurd [9], Johansen and Nielsen [22], Reisen et al. [42]). A core example -a general multiplicative model, or SARIMA (p, d, q)×(P, D, Q) s -was introduced in the book by Box and Jenkins [5]. It includes both integrated and seasonal factors: where Ψ(z) and Θ(z) are two polynomials of degrees of P and Q respectively which have roots outside the unit circle. The parameters d and D are allowed to be fractional. When |d + D| < 1/2 and |D| < 1/2, the process (1) is stationary and invertible. The paper by Porter-Hudak [41] illustrates an application of a seasonal ARFIMA model to the analysis of the monetary aggregates used by U.S. Federal Reserve. Another model of fractional integration is GARMA processes described by the equation (see Gray, Cheng and Woodward [16]) For the resent results dedicated to the statistical inference for seasonal long-memory sequences, we refer to the paper by Tsai, Rachinger and Lin [46], who developed methods of estimation of parameters in case of measurement errors. In their paper Baillie, Kongcharoen and Kapetanios [3] compared MLE and semiparametric estimation procedures for prediction problems based on ARFIMA models. Based on simulation study, they indicate better performance of MLE predictor than the one based on the two-step local Whittle estimation. Hassler and Pohle [19] (see also Hassler [18]) assess a predictive performance of various methods of forecasting of inflation and return volatility time series and show strong evidences for models with a fractional integration component.
Another type of non-stationary processes are periodically correlated, or cyclostationary, processes introduced by Gladyshev [13], which belong to the class of processes with time-dependent spectrum and are widely used in signal processing and communications (see Napolitano [37] for a review of recent works on cyclostationarity and its applications). Periodic time series are also considered as extension of seasonal models [2,4,28,38].
One of the fields of interests related to time series analysis is optimal filtering. It aims to remove the unobserved components, such as trends, seasonality or noise signal, from the observed data [1,6].
Methods of parameters estimations and filtering usually do not take into account the issues arising from real data, namely, the presence of outliers, measurement errors, incomplete information about the spectral, or model, structure etc. From this point of view, we see an increasing interest to robust methods of estimation that are reasonable in such cases (see Reisen, et al. [43], Solci at al. [45] for the examples of robust estimates of SARIMA and PAR models). Grenander [15], Hosoya [20], Kassam [24], Franke [10], Vastola and Poor [47], Moklyachuk [32,33] studied minimax extrapolation, interpolation and filtering problems for stationary sequences and processes.
This article is dedicated to the robust filtering of stochastic sequences with periodically stationary long memory multiple seasonal increments, or sequences with periodically stationary general multiplicative (GM) increments, introduced by Luz and Moklyachuk [31]. In the recent years, models with multiple seasonal and periodic patterns are in scope of interest, see Dudek [8], Gould et al. [14] and Reisen et al. [42], Hurd and Piparas [21]. This motivates us to study the time series combining a periodic structure of the covariation function as well as multiple seasonal factors. The discussed problem is a natural continuation of the researches on minimax filtering of stationary vector-valued processes, periodically correlated processes and processes with stationary increments have been performed by Moklyachuk and Masyutka [34], Moklyachuk and Golichenko (Dubovetska) [7], Luz and Moklyachuk [29,30] respectively. We also mention the works [35,36] by Moklyachuk and Sidei, who derive minimax estimates of stationary processes from observations with missed values, and the work [27] by Moklyachuk and Kozak, who have studied interpolation problem for stochastic sequences with periodically stationary increments.
The article is organized as follows. In Section 2, we recall definitions of generalized multiple (GM) increment sequence χ (d) µ,s ( ξ(m)) and stochastic sequences ξ(m) with periodically stationary (periodically correlated, cyclostationary) GM increments. The spectral theory of vector-valued GM increment sequences is discussed. Section 3 deals with the classical filtering problem for linear functionals Aξ and A N ξ which are constructed from unobserved values of the sequence ξ(m) when the spectral densities of the sequence ξ(m) and a noise sequence η(m) are known. Estimates are obtained by applying the Hilbert space projection technique to the vector sequence ξ(m) + η(m) with stationary GM increments under the stationary noise sequence η(m) uncorrelated with ξ(m). The case of non-stationary fractional integration is discussed as well. Section 4 is dedicated to the minimax (robust) estimates in cases, where spectral densities of sequences are not exactly known while some sets of admissible spectral densities are specified. We illustrate the proposed technique on the particular types of the sets, which are generalizations of the sets of admissible spectral densities described in a survey article by Kassam and Poor [25] for stationary stochastic processes.
2 Stochastic sequences with periodically stationary generalized multiple increments

Preliminary notations and definitions
Consider a stochastic sequence ξ(m), m ∈ Z, and a backward shift operator B µ with the Within the article, δ lp denotes Kronecker symbols, n l = n! l!(n−l)! .
Definition 2.1. For a stochastic sequence ξ(m), m ∈ Z, the sequence is called stochastic generalized multiple (GM) increment sequence of differentiation order d with a fixed seasonal vector s ∈ (N * ) r and a varying step µ ∈ (N * ) r or ∈ (Z \ N) r .
The multiplicative increment operator χ µ,s (B) admits the representation The explicit representation of the coefficients e γ (k) is given in [31].

Definition and spectral representation of periodically stationary GM increment
In this section, we present definition, justification and a brief review of the spectral theory of stochastic sequences with periodically stationary multiple seasonal increments.
Definition 2.3. A stochastic sequence ξ(m), m ∈ Z is called stochastic sequence with periodically stationary (periodically correlated) GM increments with period T if the mathematical expectations exist for every m, k, µ 1 , µ 2 and T > 0 is the least integer for which these equalities hold.
Θ m (k)z k and put Φ m (k) = 0 for k > q 1 , Θ m (k) = 0 for k > q 2 . Then the increment sequence is periodically stationary and allows a stationary vector representation The following theorem describes the spectral structure of the GM increment [23], [31]. where c is a vector, F (λ) is the matrix-valued spectral function of the stationary stochastic sequence χ The vector c and the matrix-valued function F (λ) are determined uniquely by the GM increment sequence χ µ,s ( ξ(m)) admits the spectral representation where is a (vector-valued) stochastic process with uncorrelated increments on [−π, π) connected with the spectral function F (λ) by the relation −π ≤ λ 1 < λ 2 < π, p, q = 1, 2, . . . , T.
Consider another vector stochastic sequence with the stationary GM increments ζ(m) = ξ(m)+ η(m), where η(m) is a vector stationary stochastic sequence, uncorrelated with ξ(m), with a spectral representation , is a stochastic process with uncorrelated increments, that corresponds to the spectral function G(λ) [17]. The stochastic stationary GM increment χ (d) µ,s ( ζ(m)) allows the spectral representation . Therefore, in the case where the spectral functions F (λ) and G(λ) have the spectral densities f (λ) and g(λ), the spectral density

GM fractional increments
In the previous subsection, we describe the GM increment sequence χ Here we consider the case of possibly fractional increment orders d i .
Within the subsection, we put the step µ = (1, 1, . . . , 1). Represent the increment operator χ where (1 − B) R0+D0 is an integrating component, R j , j = 0, 1, . . . , r, are non-negative integer numbers, 1 < s 1 < . . . < s r . Below we describe a representations d j = R j + D j , j = 0, 1, . . . , r, of the increment orders d j by stating conditions on the fractional parts D j , such that the increment sequence y(m) : is a stationary fractionally integrated seasonal stochastic sequence. For example, in case of single increment is called a fractional multiple (FM) increment sequence. Consider a generating function of the Gegenbauer polinomial: The following lemma and theorem hold true [31].
3 Hilbert space projection method of filtering
We will suppose that the following conditions are satisfied: • the minimality condition on the spectral densities f (λ) and g(λ) The latter condition is the necessary and sufficient one under which the mean square errors of estimates of functional A ξ is not equal to 0. We apply the Hilbert space estimation technique proposed by Kolmogorov [26] which can be described as a 3-stage procedure: (i) define a target element (to be estimated) of the space H = L 2 (Ω, F , P) of random variables γ which have zero mean values and finite variances, Eγ = 0, E|γ| 2 < ∞, endowed with the inner product γ 1 , γ 2 = Eγ 1 γ 2 , (ii) define a subspace of H generated by observations, (iii) find an estimate of the target element as an orthogonal projection on the defined subspace.
Note, that under conditions (12), the functional A η belongs to the space H. To find the estimate of the functional A ξ, it is sufficient to find the estimate of the functional A η: Since the functional A ζ depends on values of the stochastic sequence ζ(k) which is observed, the mean square errors of the optimal estimates A ξ and A η are connected by relation Relation (14) implies that any linear estimate A ξ of the functional A ξ allows the representation where h µ (λ) = {h p (λ)} T p=1 is the spectral characteristic of the estimate A ξ. Stage (ii). The set of the observed GM increments {χ µ,s ( η(k)) : k ≤ 0}, µ > 0 generates a closed linear subspace H 0 (ξ (d) +η (d) ) of the Hilbert space H = L 2 (Ω, F, P). The functions generate a closed linear subspace L 0 2 (p) of the Hilbert space L 2 (p). Then the relation implies a map between elements Stage (iii). The estimate A η is found as an orthogonal projection of the functional A η from the space H = L 2 (Ω, F, P) on the subspace H 0 (ξ (d) + η (d) ). This projection is characterized by the following conditions: Define the following matrix-valued Fourier coefficients: Define the vectors where the coefficients a µ (k) = a −µ (k − n(γ)), k ≥ 0, Denote by P µ , S µ and Q matrices with the matrix entries (P µ ) l,k = P µ l,k , (S µ ) l,k = S µ l+1,k−n(γ) , (Q) l,k = Q l,k , l, k ≥ 0. (12) and (13) is calculated by formula (15). The spectral characteristic h µ (λ) and the value of the mean square error ∆(f, g; A ξ) are calculated by the formulas

Theorem 3.1. A solution A ξ to the filtering problem for a vector-valued stochastic sequence ξ(m) with stationary GM increments under conditions
Proof. See Appendix.
The filtering problem for the functional A N ξ is solved directly by Theorem 3.1 by putting a(k) = 0 for k > N . To solve the filtering problem for the pth coordinate of the single vector ξ(−N ), we put a(N ) = δ p , a(k) = 0 for k = N .

Filtering of periodically stationary GM increment
Consider the filtering problem for the functionals In the same way, the functional Aη is represented as Making use of the introduced notations and statements of Theorem 3.1 we can claim that the following theorem holds true. Theorem 3.2. Let a stochastic sequence ξ(m) with periodically stationary GM increments and a stochastic periodically stationary sequence η(m) generate by formulas (27) and (29) vector-valued stochastic sequences ξ(m) and η(m) with the spectral densities matrices f (λ) = {f ij (λ)} T i,j=1 and g(λ) = {g ij (λ)} T i,j=1 . A solution Aξ to the filtering problem for the functional Aξ = A ξ under conditions (12) and (13) is calculated by formula (15) for the coefficients a(m), m ≥ 0, defined in (28). The spectral characteristic h µ (λ) = {h p (λ)} T p=1 and the value of the mean square error ∆(f ; Aξ) of the estimate Aξ are calculated by formulas (17) and (18) respectively.

The functional A M ξ can be represented in the form
, the sequence ξ(m) is determined by formula (27),   (27) and (29) vector-valued stochastic sequences ξ(m) and η(m). A solution A M ξ to the filtering problem for the functional A M ξ = A N ξ under condition (13) is calculated by formula (19) for the coefficients a(m), 0 ≤ m ≤ N , defined in (30). The spectral characteristic and the value of the mean square error of the estimate A M ξ are calculated by formulas (20) and (22) respectively.
Corollary 3.4. Let a stochastic sequence ξ(m) with periodically stationary GM increments and a stochastic periodically stationary sequence η(m) generate by formulas (27) (24) and (25) respectively. ξ(m) + η(m) at points m = 0, −1, −2, . . . are proposed in Theorem 3.1 and Corollary 3.1 in the case where the spectral density matrices f (λ) and g(λ) of the target sequence and the noise are exactly known.
In this section, we study the case where the complete information about the spectral density matrices is not available while the set of admissible spectral densities D = D f × D g is known. The minimax approach of estimation of the functionals from unobserved values of stochastic sequences is considered, which consists in finding an estimate that minimizes the maximal values of the mean square errors for all spectral densities from a class D simultaneously. This method will be applied for the concrete classes of spectral densities.
The proceed with the stated problem, we recall the following definitions [33].
Definition 4.1. For a given class of spectral densities D = D f × D g , the spectral densities f 0 (λ) ∈ D f , g 0 (λ) ∈ D g are called the least favourable densities in the class D for optimal linear filtering of the functional Aξ if the following relation holds true Definition 4.2. For a given class of spectral densities D = D f × D g the spectral characteristic h 0 (λ) of the optimal estimate of the functional Aξ is called minimax (robust) if the following relations hold true Taking into account the introduced definitions and the relations derived in the previous sections we can verify that the following lemma holds true.
Lemma 4.1. The spectral densities f 0 ∈ D f , g 0 ∈ D g which satisfy the minimality condition (13) are the least favourable spectral densities in the class D for the optimal linear filtering of the functional Aξ if the matrices P 0 µ , S 0 µ , Q 0 defined by the Fourier coefficients of the functions determine a solution of the constrained optimization problem The minimax spectral characteristic The more detailed analysis of properties of the least favorable spectral densities and the minimax-robust spectral characteristics shows that the minimax spectral characteristic h 0 and the least favourable spectral densities f 0 and g 0 form a saddle point of the function ∆(h; f, g) on the set H D × D. The saddle point inequalities where the functional ∆(h µ (f 0 , g 0 ); f, g) is calculated by the formula with The constrained optimization problem (32) is equivalent to the unconstrained optimization problem where δ(f, g|D) is the indicator function of the set D, namely δ(f |D) = 0 if f ∈ D and δ(f |D) = +∞ if (f, g) / ∈ D. The condition 0 ∈ ∂∆ D (f 0 , g 0 ) characterizes a solution (f 0 , g 0 ) of the stated unconstrained optimization problem. This condition is the necessary and sufficient condition that the point (f 0 , g 0 ) belongs to the set of minimums of the convex functional ∆ D (f, g) [10,33,34,44]. Thus, we it allows us to find the equalities for the least favourable spectral densities in some special classes of spectral densities D.
The form of the functional ∆(h µ (f 0 , g 0 ); f, g) is suitable for application the Lagrange method of indefinite multipliers to the constrained optimization problem (32). Thus, the complexity of the problem is reduced to finding the subdifferential of the indicator function of the set of admissible spectral densities. We illustrate the solving of the problem (34) for concrete sets admissible spectral densities in the following subsections. A semi-uncertain filtering problem, when the spectral density f (λ) is known and the spectral density g(λ) belongs to in class D g , is considered as well.