Minimax-robust forecasting of sequences with periodically stationary long memory multiple seasonal increments

We introduce stochastic sequences $\zeta(k)$ with periodically stationary generalized multiple increments of fractional order which combines cyclostationary, multi-seasonal, integrated and fractionally integrated patterns. We solve the problem of optimal estimation of linear functionals constructed from unobserved values of stochastic sequences $\zeta(k)$ based on their observations at points $ k<0$. For sequences with known spectral densities, we obtain formulas for calculating values of the mean square errors and the spectral characteristics of the optimal estimates of functionals. Formulas that determine the least favorable spectral densities and minimax (robust) spectral characteristics of the optimal linear estimates of 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
A variety of non-stationary and long memory time series models are introduced and investigated by researchers in the past decade (see, for example, papers by Dudek and Hurd [9], Johansen and Nielsen [24], Reisen et al. [44]). Such models are used when analyzing data which arise in different field of economics, finance, climatology, air pollution, signal processing.
Since the book by Box and Jenkins (1970), autoregressive moving average (ARMA) models integrated of order d are standard models used for time series analysis. These models are described by the equation where ε t , t ∈ Z, is a sequence of zero mean i.i.d. random variables, ψ(z), θ(z) are polynomials of p and q degrees respectively with roots outside the unit circle. This integrated ARIMA model is generalized by adding a seasonal component. A new model is described by the equation (see new edition of the book by Box and Jenkins [6] for detailes) where Ψ(z) and Θ(z) are polynomials of degrees of P and Q respectively which have roots outside the unit circle.
When an ARIMA sequence determined by equation (1) is inserted in relation (2) instead of ε t we have general multiplicative model with parameters (p, d, q) × (P, D, Q) s , d, D ∈ N * , called SARIMA (p, d, q) × (P, D, Q) s model. A good performance is shown by models which include fractional integration, that is when parameters d and D are fractional. When |d + D| < 1/2 and |D| < 1/2, a process described by equation (3) is stationary and invertible. We refer to the paper by Porter-Hudak [43] who studied seasonal ARFIMA models and applied them to the monetary aggregates used by U.S. Federal Reserve.
Closely related to fractionally integrated ARMA and GARMA processes described by equation is SARFIMA process. These processes were introduced and studied by Granger and Joyeux [15], Hosking [21], Andel [1], Gray et al. [17] in order to model longmemory stationary time series. Fractionally integrated models are a powerful tool for studying a variety of real world processes. For the resent works dedicated to statistical inference for seasonal long-memory sequences, we refer to Arteche and Robinson [2], who applied the log-periodogram and Gaussian or Whittle methods of memory parameters estimation for seasonal/cyclical asymmetric long memory processes with application to UK inflation data, and also Tsai, Rachinger and Lin [48], who developed methods of estimation of parameters in case of measurement errors. Baillie, Kongcharoen and Kapetanios [4] 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 two-step local Whittle estimation. Hassler and Pohle [20] (see also Hassler [19]) 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-stationarity is described by stochastic processes with time-dependent spectrum. A wide class of processes with time-dependent spectrum is formed by periodically correlated, or cyclostationary, processes introduced by Gladyshev [13]. These processes are widely used in signal processing and communications (see Napolitano [39] for a review of recent works on cyclostationarity and its applications). Periodic time series may be considered as an extension of a SARIMA model (see Lund [30] for a test assessing if a PARMA model is preferable to a SARMA one) and are suitable for forecasting stream flows with quarterly, monthly or weekly cycles (see Osborn [40]).
The models mentioned above are used in estimation of model's parameters and forecast issues. Meanwhile a direct application of the developed results to real data may lead to significant increasing of errors of estimates due to 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. For example, Reisen, et al. [45] proposed a semiparametric robust estimator for fractional parameters in the SARFIMA model and illustrated its application to forecast of sulfur dioxide SO 2 pollutant concentrations. Solci at al. [47] proposed robust estimates of periodic autoregressive (PAR) model.
Robust approaches are successfully applied to the problem of estimation of linear functionals from unobserved values of stochastic processes. The paper by Grenander [16] should be marked as the first one where the minimax extrapolation problem for stationary processes was formulated as a game of two players and solved. Hosoya [22], Kassam [25], Franke [10], Vastola and Poor [49], Moklyachuk [33,34] studied minimax extrapolation (forecasting), interpolation (missing values estimation) and filtering (smoothing) problems for stationary sequences and processes. Recent results of minimax extrapolation problems for stationary vector-valued processes and periodically correlated processes belong to Moklyachuk and Masyutka [35,36] and Moklyachuk and Golichenko (Dubovetska) [7] respectively. Processes with stationary increments are investigated by Moklyachuk and Luz [31,32]. We also mention works by Moklyachuk and Sidei [37,38], who derive minimax estimates of stationary processes from observations with missed values. Moklyachuk and Kozak [29] studied interpolation problem for stochastic sequences with periodically stationary increments.
In this article, we present results of investigation of stochastic sequences with periodically stationary long memory multiple seasonal increments motivated by articles by Dudek [8], Gould et al. [14] and Reisen et al. [44], who considered models with multiple seasonal patterns for inference and forecasting, and Hurd and Piparas [23], who introduced two models of periodic autoregressive time series with multiple periodic coefficients.
In Section 2, we give definition of generalized multiple (GM) increment sequence χ (d) µ,s ( ξ(m)) and introduce stochastic sequences ζ(m) with periodically stationary (periodically correlated, cyclostationary) GM increments. Such kind of non-stationary stochastic sequence combines periodic structure of the covariation function of the sequences as well as multiple seasonal factors, including the integrating one. The section also contains a short review of the spectral theory of vector-valued GM increment sequences. Section 3 deals with the classical estimation problem for linear functionals Aζ and A N ζ which are constructed from unobserved values of the sequence ζ(m) when the spectral structure of the sequence ζ(m) is known. Estimates are obtained by representing the sequence ζ(m) as a vector sequence ξ(m) with stationary GM increments and applying the Hilbert space projection technique. An approach to forecasting in the presence of non-stationary fractional integration is discussed in Section 4. Section 5 contains examples of forecasting of particular models of time series. In Section 6, we derive the minimax (robust) estimates in cases, where spectral densities of sequences are not exactly known while some sets of admissible spectral densities are specified which are generalizations of the corresponding sets of admissible spectral densities described in a survey article by Kassam and Poor [26] for stationary stochastic processes.
2 Stochastic sequences with periodically stationary generalized multiple increments

Definition and spectral representation of a 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. This type of stochastic sequences will allow us to deal with a wide range of non-stationarity in time series analysis. Consider a stochastic sequence {η(m), m ∈ Z}. By B µ denote a backward shift operator with the step µ ∈ Z, such that B µ η(m) = η(m − µ); B := B 1 . Recall the following definition [32,42,52].
where n l = n! l!(n−l)! , is called stochastic nth increment sequence with a step µ ∈ Z.
Note that in Definition 2.1, the step parameter µ is not fixed and varies over the set Z. The introduced increment (5) is applicable for describing the integrated stochastic sequence (1). The varying step µ provides a flexibility of the integrated processes. For instance, let a sequence x m satisfy the equation is stationary as a sum of stationary 1-step increments. To deal with seasonal time series (2) we need to extend definition of stochastic increment sequence as follows.
Definition 2.2. For a given stochastic sequence {η(m), m ∈ Z}, the sequence is called stochastic seasonal increment sequence with a fixed seasonal parameter s ∈ N * = N \ {0} and a varying step µ ∈ Z.
We mention the following properties of the seasonal increment sequence η (n) s (m, µ), which will be used for proving Theorem 2.2: where {A l , l = 0, 1, 2, . . . , (µ − 1)n} are coefficients from the representation General multiplicative model (3) [6] indicates the necessity of dealing with increments of different seasonal parameters. Moreover, for each season factor at each differencing order it is possible to make different steps by applying the operator (1 − B s µ1 ) · . . . · (1 − B s µn ) instead of (1 − B s µ ) n . Thus, the following generalization is reasonable.
where all roots of polynomials φ(z), θ(z), Φ i (z), Θ i (z) lie outside the unit circle, x m is stationary in this case and we can define a GM increment sequence χ Let γ denotes a triple (µ, s, d). For i = 1, 2, . . . , r, j ∈ Z define coefficients M j i := j µisi and I j i := I{j mod µ i s i = 0}, where I{·} is the indicator function, and notations n i : Denote the maximun of two numbers as x ∨ y and the minimum as x ∧ y.
Proof. See Appendix.
µ,s (η(m)) generated by a stochastic sequence {η(m), m ∈ Z} is wide sense stationary if the mathematical expectations exist for all m 0 , m, µ, µ 1 , µ 2 and do not depend on m 0 . The function c  µ,s (η(m)) by (9) is called stochastic sequence with stationary GM increments (or GM increment sequence of order d).
The following theorem is a generalization of the corresponding theorem for stochastic increment sequence η (d) (m, µ) [32,51].    Note that by the spectral function and the spectral density of a stochastic sequence with stationary GM increments, we will call the spectral function and the spectral density of the corresponding stationary GM increment sequence. Representation (11) and the Karhunen theorem [27,11] imply the spectral representation of the stationary GM increment sequence χ (d) µ,s (η(m)): where Z η (d) (λ) is a stochastic process with uncorrelated increments on [−π, π) connected with the spectral function F (λ) by the relation Finally, we are ready to give a definition of periodically stationary GM increment sequence.
Definition 2.5. A stochastic sequence {ζ(m), m ∈ Z} is called stochastic sequence with periodically stationary (periodically correlated) GM increments with period T if the mathematical expectations T s (m, k; µ 1 , µ 2 ) exist for every m, k, µ 1 , µ 2 and T > 0 is the least integer for which these equalities hold.
It follows from Definition 2.5 that the sequence forms a vector-valued sequence ξ(m) = {ξ p (m)} p=1,2,...,T , m ∈ Z with stationary GM increments as follows: µ,s (ξ p (m)) is the GM increment of the p-th component of the vectorvalued sequence ξ(m). The following theorem describes a spectral representation of the sequence ξ(m).
Representation (17) is called the canonical moving average representation of the stochastic stationary GM increment sequence χ The spectral function F (λ) of a stationary GM increment sequence χ (d) µ,s ( ξ(m)) which admits the canonical representation (17) has the spectral density f (λ) = {f ij (λ)} T i,j=1 admitting the canonical factorization where the function Φ(z) = ∞ k=0 ϕ(k)z k has analytic in the unit circle {z : Then the following relation holds true: We will use the one-sided moving average representation (17) and relation (19) for finding the mean square optimal estimates of unobserved values of vector-valued sequences with stationary GM increments.
3 Hilbert space projection method of forecasting

Forecasting of vector-valued stationary GM increment
Consider a vector-valued stochastic sequence with stationary GM increments ξ(m) constructed from the sequence ζ(m) with the help of transformation (13). Let the stationary GM increment sequence χ µ,s (ξ p (m))} T p=1 has an absolutely continuous spectral function F (λ) and the spectral density f (λ) = {f ij (λ)} T i,j=1 . Without loss of generality we will assume that Eχ µ,s ( ξ(m)) = 0 and µ > 0.
Consider the problem of mean square optimal linear estimation of the functionals which depend on unobserved values of the stochastic sequence ξ with stationary GM increments. Estimates are based on observations of the sequence ξ(k) at points k = −1, −2, . . .. We will suppose that the following conditions are satisfied: • the minimality condition on the spectral density f (λ) The latter one is the necessary and sufficient condition under which the mean square errors of estimates of functionals A ξ and A N ξ are not equal to 0. The classical Hilbert space estimation technique proposed by Kolmogorov [28] can be described as a 3-stage procedure: (i) define a target element of the space H = L 2 (Ω, F , P) to be estimated, (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.
Stage i. Neither functional A ξ nor A N ξ belongs to the space H. With the help of the following lemma and the corresponding corollary, we describe representations of these functionals as sums of functionals with finite second moments belonging to H and functionals depending on observed values of the sequence ξ(k) ("initial values").
Proof. See Appendix.
D µ N is the linear transformation determined by an infinite matrix with the entries So that, Lemma 3.1 provides a representation of the functional A ξ as a sum of an element Bχ ξ from the space H = L 2 (Ω, F , P) under conditions (21) - (22) and linear combination V ξ of a finite number of initial values ξ(k), k = −1, −2, . . . , −n(γ), which are observed. Thus, the following equality hold true Denote by ∆(f, A ξ) := E|A ξ − A ξ| 2 the mean square error of the optimal estimate A ξ of the functional A ξ and let ∆(f, Bχ ξ) := E|Bχ ξ − Bχ ξ| 2 denote the mean square error of the optimal estimate Bχ ξ of the functional Bχ ξ. Then Thus, we have defined the functional Bχ ξ to be used in finding the optimal linear estimate of the functional A ξ.
implies a relation between elements χ The spectral representation of the functional Bχ ξ can be written in the form Thus, at stage iii, the problem is equivalent to finding a projection of the

Relation (26) implies that every linear estimate A ξ of the functional A ξ can be written in the form
where h µ (λ) = {h p (λ)} T p=1 is the spectral characteristic of the estimate B ξ, which is a projection of the element B µ (e iλ ) This estimate is characterized by the following conditions: Condition (30) implies the following relation holding true for all k −1 Thus, the spectral characteristic of the estimate Bχ ξ can be represented in the form and c(k) = {c p (k)} T p=1 , k 0 are unknown coefficients to be found. Condition (29) implies the following representation of the spectral characteristic h µ (λ) which allows us to write the relations (32) Next we define the matrix-valued Fourier coefficients and rewrite relation (32) as a system of linear vector equations determining the unknown coefficients c µ (k), k ≥ 0. This system can be presented in the matrix form where To show that operator F µ is invertible we note that the problem of projection of the element B ξ of the Hilbert space H on the closed convex set H 0− ( ξ (d) µ ) has a unique solution for each non-zero coefficients { a(0), a(1)), a(2), . . .}, satisfying conditions (21) - (22). Therefore, equation (34) has a unique solution for each vector D µ a, which implies existence of the inverse operator F −1 µ . Therefore, coefficients c µ (k), k ≥ 0, which determine the spectral characteristic h µ (λ), can be calculated as where The spectral characteristic h µ (λ) of the estimate Bχ ξ is calculated by the formula (36) The value of the mean square error of the estimate A ξ is calculated by the formula Next consider the problem in the case where the GM incremental sequence of the stochastic sequence ξ(m) admits moving-average representation (17) and its spectral density f (λ) = {f ij (λ)} T i,j=1 admits the canonical factorization (18), (19), namely where Formulas for calculating the spectral characteristic h µ (λ) and the value of the mean square error ∆(f ; A ξ) of the estimate A ξ can be presented in terms of the function Ψ µ (e −iλ ) and the factorization coefficients ϕ µ (k), k = 0, 1, 2, . . . . One can directly check that conditions (29) and (30) are satisfied by the function where and A is a linear symmetric operator which is determined by a matrix with the entries A(k, j) = a(k + j), k, j ≥ 0. The defined operators D µ A and A are compact under conditions (21) - (22). Then the value of the mean square error is calculated by the formula The derived results are summarized in the following theorem.
which satisfies the minimality condition (23). Let coefficients a(j), j 0 satisfy conditions (21) - (22). Then the optimal linear estimate A ξ of the functional A ξ based on observations of the sequence ξ(m) at points m = −1, −2, . . . is calculated by formula (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 (36) and (37) respectively.
In the case where the spectral density f (λ) admits the canonical factorization (38) the spectral characteristic and the value of the mean square error of the optimal estimate Aξ can be calculated by formulas (39) and (40) respectively.

Estimates of functional
is calculated by the formula where and D µ N is defined in Corollary 3.1. The value of the mean square error of the estimate A N ξ is In the case where the spectral density f (λ) admits the canonical factorization (38) the spectral characteristic can be calculated as where (1), . . . , ϕ µ (N )); A N is a linear operator determined by the coefficients a(k), k = 0, 1, . . . , N , as follows: The value of the mean square error is calculated by the formula Thus, the following theorem holds true.
, m ∈ Z} be a stochastic sequence which determine a stationary stochastic GM increment sequence χ (n) µ,s ( ξ(m)) with the spectral density matrix f (λ) which satisfies the minimality condition (23). The optimal linear estimate A N ξ of the functional A N ξ based on observations of the sequence ξ(m) at points m = −1, −2, . . . is calculated by formula (41). The spectral characteristic h µ,N (λ) = {h µ,N,p (λ)} T p=1 and the value of the mean square error ∆(f ; A N ξ) are calculated by formulas (42) and (43) respectively. In the case where the spectral density f (λ) admits the canonical factorization (38) the spectral characteristic h µ,N (λ) and the value of the mean square error of the estimate A N ξ can be calculated by formulas (44) and (45) respectively.
For the problem of the mean square optimal estimate of the unobserved value A N,p ξ = ξ p (N ) = ξ(N )δ p , p = 1, 2, . . . , T , N ≥ 0 of the stochastic sequence ξ(m) with GM stationary increments based on its observations at points m = −1, −2, . . . we have the following corollary from Theorem 3.2.
Corollary 3.2. The optimal linear estimate ξ p (N ) of the value ξ p (N ), p = 1, 2, . . . , T , N ≥ 0, of the stochastic sequence with GM stationary increments from observations of the sequence ξ(m) at points m = −1, −2, . . . is calculated by formula The spectral characteristic h µ,N,p (λ) of the estimate is calculated by the formula The value of the mean square error of the estimate ξ p (N ) is calculated by the formula In the case where the spectral density f (λ) admits canonical factorization (38), and the condition min i=1,r µ i s i > N is satisfied, the spectral characteristic and the value of the mean square error of the estimate ξ p (N ) can be calculated by the formulas and Remark 3.1. Since for all d ≥ 1 and µ ≥ 1 the condition which can be calculated by the formula For this reason the following relation holds true:

Forecasting of periodically stationary GM increment
Consider the problem of mean square optimal linear estimation of the functionals (56) Making use of the introduced notations and statements of Theorem 3.1 we can claim that the following theorem holds true.  (36) and (37) respectively. In the case where the spectral density matrix f (λ) admits the canonical factorization (38), the spectral characteristic and the value of the mean square error of the estimate Aξ can be calculated by formulas (39) and (40) respectively.

The functional A M ζ can be represented in the form
, the sequence ξ(m) is determined by formula (55), Making use of the introduced notations and statements of Theorem 3.2 we can claim that the following theorem holds true.  (42) and (43) respectively. In the case where the spectral density matrix f (λ) admits the canonical factorization (38), then the spectral characteristic h µ,N (λ) and the value of the mean square error of the estimate A M ζ can be calculated by formulas (44)

Forecasting of GM fractional increments
In the previous section, we solved the forecasting problem for the increment sequence χ where (1 − B) R0+D0 is the integrating component, R j , j = 0, 1, . . . , r, are nonnegative integer numbers, 1 < s 1 < . . . < s r . The goal is to find representations d j = R j + D j , j = 0, 1, . . . , r, of the increment orders under some conditions on the fractional parts D j , such that the increment sequence y(m) : to be a stationary fractionally integrated seasonal stochastic sequence. For example, in case of single increment pattern (1 − B s * ) R * +D * this condition is |D * | < 1/2.
We will call a sequence χ   Lemma 4.2 shows that a multiple seasonal increment sequence can be represented as the following k-factor Gegenbauer sequence In the case where ξ(m) is an ARM A(p, q) sequence, the model x(m) defined by (59) is called k-factor GARM A(p, d i , u i , q) sequence. It is stationary and invertible if |d i | < 1/2 for |u i | < 1 and |d i | < 1/4 for |u i | = 1. If additionally d i > 0, then the model exhibits a long memory behavior [50]. The function (1 − 2u i B + B 2 ) −di is a generating function of the Gegenbauer polynomial: Thus, denoting k * = |M|, we obtain G − k * (m) = 0≤n1,...,n k * ≤m,n1+...+n k * =m ν∈M The derived representations of the increment operator χ (D) s (B) imply the following theorem.
In the following remarks we provide some additional details with the help of which we can use theorems proposed in the previous section in finding solution of the forecasting problem for stochastic sequences with periodically stationary (periodically correlated) FM increments.

It is stationary under conditions
3, D 1 = 0.2 and the process exhibits a long-memory behavior.

and admits a representation
We find the spectral characteristic h (1,1),2 (λ) of the estimate A 2 ξ using Theorem 3.2 as well as remarks to Theorem 4.1. First we obtain The optimal estimate of the functional A 2s ξ is calculated by the formula The value of the mean square error of the estimate is calculated by the formula In a particular case d 0 = 1, d 1 = 1 and, respectively, D 0 = 0, D 1 = 0, we have χ (D) (1,u) (e −iλ ) ≡ 1, and G ± k * (k) = 0 for k ≥ 1. In this case the estimate of the functional A 2s ξ and the value of the its mean square error are calculated by the formulas

Minimax (robust) method of forecasting
Values of the mean square errors and spectral characteristics of the optimal estimates of functionals A ξ and A N ξ constructed from unobserved values of stochastic sequence ξ(m) which determine a stationary stochastic GM increment sequence χ µ,s ( ξ(m)) with the spectral density matrix f (λ) based on its observations ξ(m) at points m = −1, −2, . . . can be calculated by formulas (37), (36) and (43), (42) respectively, in the case where the spectral density matrix f (λ) is exactly known. In the case where the spectral density f (λ) admits the canonical factorization (38), formulas (71), (39) and (45), (44) are derived for calculating values of the mean square errors and spectral characteristics, respectively.
In many practical cases, however, complete information about the spectral density matrix is impossible while some sets D of admissible spectral densities can be defined. In this case the minimax method of estimation of functionals from unobserved values of stochastic sequences is reasonable. This method consists in finding an estimate that minimizes the maximal values of the mean square errors for all spectral densities from a given class D of admissible spectral densities simultaneously. Definition 6.1. For a given class of spectral densities D a spectral density f 0 (λ) ∈ D is called the least favourable in D for the optimal linear estimation of the functional A ξ if the following relation holds true: Definition 6.2. For a given class of spectral densities D a spectral characteristic h 0 (λ) of the optimal linear estimate of the functional Aξ is called minimax-robust if the following conditions are satisfied Taking into account the introduced definitions and the relations derived in the previous sections we can verify that the following lemmas hold true.
determines a solution to the constrained optimization problem The minimax spectral characteristic h 0 = h µ (f 0 ) is calculated by formula (36) if h µ (f 0 ) ∈ H D .

Lemma 6.2.
A spectral density f 0 (λ) ∈ D which admits the canonical factorization (38) is the least favourable density in the class D for the optimal linear estimation of the functional A ξ based on observations of the sequence ξ(m) at points m = −1, −2, . . . if coefficients {ϕ 0 (k) : k ≥ 0} of the canonical factorization of the spectral density f 0 (λ) determine a solution to the constrained optimization problem The minimax spectral characteristic h 0 = h µ (f 0 ) is calculated by formula (39) if h µ (f 0 ) ∈ H D .
of the spectral density f 0 (λ) determine a solution to the constrained optimization problem For more detailed analysis of properties of the least favorable spectral densities and the minimax-robust spectral characteristics we observe that the minimax spectral characteristic h 0 and the least favourable spectral density f 0 form a saddle point of the function ∆(h; f ) on the set H D × D. The saddle point inequalities ∈ H D and f 0 is a solution of the constrained optimization problem where the functional ∆(h µ (f 0 ); f ) is calculated by the formula or by the formula in the case where the spectral density admits the canonical factorization (38). The constrained optimization problem (69) is equivalent to the unconstrained optimization problem where δ(f |D) is the indicator function of the set D, namely δ(f |D) = 0 if f ∈ D and δ(f |D) = +∞ if f / ∈ D. A solution f 0 of this unconstrained optimization problem is characterized by the condition 0 ∈ ∂∆ D (f 0 ), which is the necessary and sufficient condition under which a point f 0 belongs to the set of minimums of the convex functional ∆ D (f ) [10,35,33,46]. This condition makes it possible to find the least favourable spectral densities in some special classes of spectral densities D.
The form of the functional ∆(f ) allows us to apply the Lagrange method of indefinite multipliers for investigating the constrained optimization problem (69). The complexity of optimization problem is determined by the complexity of calculating the subdifferentials of the indicator functions of sets of admissible spectral densities.

Least favorable spectral density in classes D 0
Consider the forecasting problem for the functional A ξ which depends on unobserved values of a sequence ξ(m) with stationary GM increments based on observations of the sequence at points m = −1, −2, . . . under the condition that sets of admissible spectral densities D k 0 , k = 1, 2, 3, 4 are defined as follows: where p, p k , k = 1, T are given numbers, P, B 1 , are given positive-definite Hermitian matrices.
From the condition 0 ∈ ∂∆ D (f 0 ) we find the following equations which determine the least favourable spectral densities for these given sets of admissible spectral densities.
For the first set of admissible spectral densities D 1 0 we have equation where α is a vector of Lagrange multipliers.
For the second set of admissible spectral densities D 2 0 we have equation where α 2 is a Lagrange multiplier.
For the third set of admissible spectral densities D 3 0 we have equation where α 2 k are Lagrange multipliers, δ kl are Kronecker symbols. For the fourth set of admissible spectral densities D 4 0 we have equation where α 2 is a Lagrange multiplier.
In the case where the spectral density admits the canonical factorization (38) we have the following equations, correspondingly The following theorem holds true.
Here the spectral densities V (λ), U (λ) are known and fixed, q, q k , k = 1, T are given numbers, Q, B 2 are given positive definite Hermitian matrices. From the condition 0 ∈ ∂∆ D (f 0 ) we find the following equations which determine the least favourable spectral densities for these given sets of admissible spectral densities.
For the second set of admissible spectral densities D U For the third set of admissible spectral densities D U where β 2 k are Lagrange multipliers, δ kl are Kronecker symbols, γ 1k (λ) ≤ 0 and In the case where the spectral density admits the canonical factorization (38) we have the following equations, correspondingly , (87) The following theorem holds true.

Least favorable spectral density in classes D ε
Consider the forecasting problem for the functional A ξ which depends on unobserved values of a sequence ξ(m) with stationary GM increments based on observations of the sequence at points m = −1, −2, . . . under the condition that sets of admissible spectral densities D k ε , k = 1, 2, 3, 4 are defined as follows: Here f 1 (λ) is a fixed spectral density, W (λ) is an unknown spectral density, p, p k , k = 1, T , are given numbers, P is a given positive-definite Hermitian matrices.
From the condition 0 ∈ ∂∆ D (f 0 ) we find the following equations which determine the least favourable spectral densities for these given sets of admissible spectral densities.
For the first set of admissible spectral densities D 1 ε we have equation where α 2 is Lagrange multiplier, γ 1 (λ) ≤ 0 and γ 1 For the second set of admissible spectral densities D 2 ε we have equation For the third set of admissible spectral densities D 3 ε we have equation For the fourth set of admissible spectral densities D 4 ε we have equation where α is a vector of Lagrange multipliers, Γ(λ) ≤ 0 and Γ(λ In the case where the spectral density admits the canonical factorization (38) we have the following equations, correspondingly The following theorem holds true.  (38), respectively), the constrained optimization problem (64) and restrictions on densities from the corresponding classes D k ε , k = 1, 2, 3, 4. The minimax-robust spectral characteristic of the optimal estimate of the functional A ξ is determined by the formula (36).

Least favorable spectral density in classes D 1δ
Consider the forecasting problem for the functional A ξ which depends on unobserved values of a sequence ξ(m) with stationary GM increments based on observations of the sequence at points m = −1, −2, . . . under the condition that the sets of admissible spectral densities D k 1δ , k = 1, 2, 3, 4 are defined as follows: Here f 1 (λ) is a fixed spectral density, δ, δ k , k = 1, T , δ j i , i, j = 1, T , are given numbers.
From the condition 0 ∈ ∂∆ D (f 0 ) we find the following equations which determine the least favourable spectral densities for these given sets of admissible spectral densities.
For the first set of admissible spectral densities D 1 1δ we have equation where β 2 is Lagrange multiplier, |γ 2 (λ)| ≤ 1 and For the second set of admissible spectral densities D 2 1δ we have equation where β 2 k are Lagrange multipliers, γ 2 k (λ) ≤ 1 and For the third set of admissible spectral densities D 3 1δ we have equation where β 2 is a Lagrange multiplier, |γ ′ 2 (λ)| ≤ 1 and For the fourth set of admissible spectral densities D 4 1δ we have equation where β ij are Lagrange multipliers, |γ ij (λ)| ≤ 1 and In the case where the spectral density admits the canonical factorization (38) we have the following equations, correspondingly The following theorem holds true.  (38), respectively), the constrained optimization problem (64) and restrictions on densities from the corresponding classes D k 1δ , k = 1, 2, 3, 4. The minimax-robust spectral characteristic of the optimal estimate of the functional A ξ is determined by the formula (36).

Conclusions
In this article, we present results of investigation of stochastic sequences with periodically stationary long memory multiple seasonal increments. We give definition of generalized multiple increment sequence and introduce stochastic sequences ζ(m) with periodically stationary (periodically correlated, cyclostationary) generalized multiple increments. These non-stationary stochastic sequences combine periodic structure of covariation functions of sequences as well as multiple seasonal factors, including the integrating one. A short review of the spectral theory of vector-valued generalized multiple increment sequences is presented. We describe methods of solution of the classical forecasting problem for linear functionals which are constructed from unobserved values of a sequence with periodically stationary generalized multiple increments in the case where the spectral structure of the sequence is exactly known. Estimates are obtained by representing the sequence under investigation as a vectorvalued sequence with stationary generalized multiple increments and applying the Hilbert space projection technique. An approach to forecasting in the presence of non-stationary fractional integration is discussed. Examples of solution of the forecasting problem for particular models of time series are proposed. The minimax-robust approach to forecasting problem is applied in the case of spectral uncertainty where densities of sequences are not exactly known while, instead, sets of admissible spectral densities are specified. We propose a representation of the mean square error in the form of a linear functional in L 1 with respect to spectral densities, which allows us to solve the corresponding constrained optimization problem and describe the minimax (robust) estimates of the functionals. Relations are described which determine the least favourable spectral densities and the minimax spectral characteristics of the optimal estimates of linear functionals for a collection of specific classes of admissible spectral densities.

Proof of Theorem 2.2
We follow the idea proposed by Yaglom [51] for continuous time stationary increments. Consider a GM increment sequence with one seasonal factor χ The latter equality implies where the function β (d) (iu) is to be chosen in the way that the integrals are defined for λ ∈ [−π, π) and converge at the neighborhoods of the points cos su = 1, u ∈ [−π, π], namely, the points u = 2πk/s for |k| ≤ s/2, k ∈ Z. Thus, we choose the function β (d) (iu) = The function F (λ) is real-valued non-decreasing bounded function defined on [−π, π), such that F (0) = 0. Consider the function F (λ) as left-continuous. This function is the one stated in the theorem.
Proof of Lemma 3.1 Using Definition 2.3 we obtain the formal equality and