Spectral Density Estimation for Nonstationary Data With Nonzero Mean Function

Abstract We introduce a new approach for nonparametric spectral density estimation based on the subsampling technique, which we apply to the important class of nonstationary time series. These are almost periodically correlated sequences. In contrary to existing methods, our technique does not require demeaning of the data. On the simulated data examples, we compare our estimator of spectral density function with the classical one. Additionally, we propose a modified estimator, which allows to reduce the leakage effect. Moreover, in the supplementary materials, we provide a simulation study and two real data economic applications. Supplementary materials for this article are available online.


Introduction
Spectral theory for stationary time series is well developed (see Priestley 1981;Žurbenko 1986;Stoica and Moses 2005, and many others).However, in the real data applications one deals often with nonstationary phenomena.In such a case, statistical tools designed for stationary processes should be appropriately modified or generalized.If this is not possible, then a new approach needs to be proposed.
In this article, we focus on the spectral density estimation problem in the nonstationary setting.Classical approach for spectral density estimation requires removing the mean function from the data.However, as we discuss below, this may be very challenging in the nonstationary case.Thus, it is crucial to propose a new method that provides consistent estimates and does not assume the nullity of the mean function.
In the sequel, we introduce and develop a novel methodology for spectral densities estimation for nonstationary mixing sequences with nonzero mean function.Our main results cover the class of almost periodically correlated (APC) time series, which is an important class of nonstationary time series containing among others periodically correlated (PC) and covariance stationary sequences.APC and PC processes are broadly used to model data with cyclic features.Many interesting examples concerning, for example, telecommunication, vibroacoustics, mechanics, economics and climatology can be found in Antoni (2009), Gardner, Napolitano, and Paura (2006), Hurd and Miamee (2007), and Napolitano (2012Napolitano ( , 2016)).To be precise, a time series {X t , t ∈ Z} with finite second moments is called PC if it has periodic mean and covariance functions, while in APC case the mean and covariance functions are almost periodic functions of the time argument.A function f (t) : Z −→ R is called almost periodic if for any > 0 there exists an integer L > 0 such that among any L consecutive integers there is an integer p such that sup t∈Z |f (t + p ) − f (t)| < (see Corduneanu 1989).
Problem of spectral density estimation for PC and APC sequences is not new (see Nematollahi and Rao 2005;Lii and Rosenblatt 2006;Hurd and Miamee 2007;Napolitano 2012Napolitano , 2016)), but as we mentioned before all methods assume that EX t ≡ 0. In the case of stationary and PC time series this is not a restrictive condition, because the data can be demeaned (see, e.g., Brillinger 2001;Hurd and Miamee 2007).However, for APC time series this problem is much more difficult.In general case, it is not easy to estimate an almost periodic mean function.Moreover, in the presence of an almost periodic mean, the spectral densities and cyclic spectra contain Dirac impulses as a consequence of the presence of an almost periodic component in the cyclic autocovariance (see Napolitano 2016).Estimation of the spectral densities or the cyclic spectrum without removing the mean, results in spikes in correspondence of the ideal Dirac deltas.In such a case, Napolitano and Spooner (2000) proposed a median filtering technique to remove the spikes.However, for now the theoretical properties of this approach were not studied in the literature.
Our idea for spectral densities estimation is based on the subsampling method, which allows to approximate the asymptotic variance of the rescaled statistics of interest.Let us indicate that the general consistency result concerning subsampling estimator of the asymptotic variance that we present does not require periodicity or almost periodicity of the data.However, we apply it to an APC process to provide at the same time an application that is interesting and important from practical point of view.For that purpose we use results from Lenart (2013).Lenart considered estimation problem for the Fourier coefficients of the mean function of an APC time series.He showed that their rescaled estimators are asymptotically normal and that elements of the asymptotic covariance matrix are linear functions of the spectral densities of the considered process.Applying subsampling to estimate elements of this matrix we get consistent estimators of the spectral densities.Moreover, our subsampling estimator of the spectral density functions turns out to be a generalized type of a periodogram-subsampling estimator and it is associated with a classical spectral density estimator with the Bartlett window.
Article is organized as follows.In Section 2, the known methods for spectral density estimation in stationary and almost periodic case are summarized.The main results are provided in Section 3. In particular, the subsampling estimator of the spectral density functions of an APC time series with the nonzero mean function is constructed and its consistency is shown.Moreover, to reduce the leakage effect the modification of this estimator is discussed.Finally, its relationship with the periodogramsubsampling estimator and the classical spectral density estimator with the Bartlett window is derived.Results of the performed simulation study, real data example and all proofs are presented in the supplementary materials.

Spectral Density Estimation for APC Time Series
Let {X t , t ∈ Z} be an APC time series.Its mean function μ(t) = E(X t ) and the autocovariance function B(t, τ ) = cov(X t , X t+τ ) (for any fixed τ ∈ Z) can be represented as Fourier series where the sets Symbol ∼ indicates some association between the Fourier series and the almost periodic function.The sup-norm used in the definition of an almost periodic function (see Section 1) results in the uniform convergence of the Fourier series (1) to μ(t) and B(t, τ ), respectively.For more details, we refer the reader to Napolitano (2012).
In the stationary case, that is, when = {0}, the consistent estimator of the spectral density function is obtained by tapering the covariance estimator or equivalently by smoothing the periodogram (see, e.g., Priestley 1981;Brillinger 2001).Let us recall that when a PC process with known period d is considered, then the spectral density function P(ν, ω) can be nonzero only on the 2d Hurd 1991;Dehay and Hurd 1994).As a result the consistent estimator of the spectral density functions can be obtained by smoothing the shifted periodogram along the support lines or equivalently by weighting the estimators of the Fourier coefficients of the autocovariance function (Hurd and Miamee 2007, chap. 10.2).On the other hand, in Nematollahi and Rao (2005) to construct an estimator of the spectral density matrix, the relation between the PC time series and the multivariate stationary time series is used.Nematollahi and Rao proposed to estimate the spectral matrix of the stationary vector time series and then to use it to construct the estimator of the spectral density matrix of the considered PC time series.Finally, in Lii and Rosenblatt (2006), the estimation of the frequencies in almost periodic autocovariance structure of a zero mean Gaussian harmonizable APC time series is considered.In Napolitano (2016), the review of known results for the APC processes is provided.Two estimators of the spectral density functions are discussed: the frequency-smoothed cyclic periodogram and the time-smoothed cyclic periodogram (see Napolitano 2016, sec.2).
Independently on the type of the process (stationary, PC, APC) to construct the aforementioned estimators of the spectral density function, the nullity of the mean function is assumed.If the considered process is not zero mean, before one estimates the spectral density, then the data should be demeaned.This can be easily achieved when the mean function is constant or periodic with the known period length.It is also possible when one deals with an APC time series and the set is finite and known.Then, having sample (X 1 , . . ., X n ), one can apply the estimator of the mean function Finally, the estimator of the spectral density P(ν, ω) can be obtained by replacing μ(t) by its estimate μ n (t) (see Lenart 2011Lenart , 2016)).According to our knowledge such estimator was not considered in the literature in APC case when the set is unknown and needs to be estimated.It is worth to note that in the general APC case, when μ(t) is an almost periodic function with the unknown set , it is very difficult to remove entirely the mean function from the data.Moreover, even in purely periodic case one may deal with so called hidden periodicities (Yavorskyj et al. 2012) and in consequence estimation of the mean function can be challenging.In Napolitano (2016, see (2.9)), it was shown that if the almost periodic mean is not removed from the data, the classical estimates of the spectral densities for zero mean case contain Dirac impulses, which can be observed as additional spikes.The shape of the spikes is related to the length of the frequency-smoothing window if one uses the frequency-smoothing cyclic periodogram method for the estimation.In such a case, the estimator of spectral density function is not consistent.To solve this problem in Napolitano and Spooner (2000), a median filtering technique is proposed, but so far the theoretical properties of this approach were not studied in the literature.In a very special case, when the considered process is a sum of periodic mean component and a stationary time series, one can estimate spectrum without removing the periodic component.The overview of possible techniques can be found in Li (2016, see chap. 7.2).However, these methods are not suitable for a general PC time series (with a periodic autocovariance function) or more complex nonstationary time series like APC processes.For instance, they assume that the spectral density function is nonzero, which in the APC case means that we know the set .This is a restrictive condition that we do not use for our results.
In the next section, we introduce a new method for the spectral density function P(ν, ω) (see ( 2)) estimation on the bifrequency square (ν, ω) ∈ [0, 2π) 2 in a case when the almost periodic mean function is unknown, that is, when set is unknown.Our approach does not require estimation of the mean function and can be applied directly to the raw data.Let us indicate that the accurate estimation of the mean function is crucial for the second-order analysis of signals, that is, analysis of the characteristics based on the second-order moments.In the case of PC time series for that purpose one needs to know the period length d.When d is unknown (i.e., when the set is unknown), it needs to be estimated.Usually, the firstorder frequencies (the true frequencies of the mean function) are much stronger than the second-order frequencies (the true frequency of the autocovariance function).This means that if λ 1 is the first-order true frequency and λ 2 the second-order true frequency, usually for the corresponding Fourier coefficients of the autocovariance function we have |a(λ 1 , τ )| > |a(λ 2 , τ )|.Of course, if the mean function is known to be periodic and the period length d is known, the simplest and the most natural approach is the standard one, that is, demeaning the data and then performing the analysis.But if d is unknown our approach allows to avoid its estimation and consequently, we are not dependent on the quality of this estimate.Additionally, our method can also be applied when the mean function is almost periodic.Moreover, in the real data applications, it may happen that process that seems to be periodic, in reality is not purely periodic.For instance, when the sample size grows, one may observe irregular cyclities (see Gardner 2019) and then modeling data with the PC process is no longer correct.An interesting example is also the electrocardiogram signal, which can be considered as PC only if its length has a few seconds (Napolitano 2017).

Main Results
In the first part of this section, we propose the subsampling covariance matrix estimator for general weakly dependent time series and show its mean square consistency.In the next step, we use this general result in the construction of the spectral density function estimator for an APC time series in the case when the sets and are unknown.
From now on by Re(z) and Im(z), we denote the real and the imaginary part of a complex number z. Symbol (•) denotes transpose of a vector.By N s (a, ) we denote the s-variate Gaussian distribution with mean a and covariance matrix .Finally, for any random variable Z and real number p > 0, Z p = E|Z| p 1/p .

General Subsampling Consistency Result
The idea of subsampling (see Politis, Romano, and Wolf 1999) is to approximate the sampling distribution of the statistics using subsamples of the data.Let {X t , t ∈ Z} be a considered time series and (X 1 , . . ., X n ) be an observed sample.Moreover, let r be any positive integer.By θ = θ (1) , θ (2) , . . ., θ (r) ∈ R r and , we denote an unknown parameter and its estimator based on the sample (X 1 , X 2 , . . ., X n ).We assume that there exists a nondegenerate limiting distribution of τ n ( θ n − θ ), where τ n is a normalizing sequence.Finally, θ n,b,t is a subsampling version of θ n based on the subsample The important feature of subsampling (and its recent generalizations) is that it provides valid quantiles for confidence intervals and hypothesis tests in stationary and nonstationary models under very mild regularity conditions (see, e.g., Politis, Romano, and Wolf 1999;Dehay, Dudek, and Leśkow 2014;Dudek and Lenart 2017;Lenart 2011Lenart , 2013Lenart , 2018;;Tewes, Politis, and Nordman 2019).Moreover, to obtain subsampling consistency one does not need to know the exact form of the limiting distribution of τ n ( θ n − θ ).However, in this article, we are only interested in a subsampling approximation of the elements of the covariance matrix of this limiting distribution.By = [σ ij ] r×r , we denote this limiting covariance matrix.
Recall that for fixed 1 ≤ i, j ≤ r and any sequence b such that b → ∞, b/n → 0 the estimator is a subsampling estimator of θ (i) , see Politis, Romano, and Wolf (1999, sec. 3.8.1 and 4.6).However, in this article, we define some more general subsampling estimator of σ ij , which is based on the generalized subsampling procedure introduced by Lenart (2018).For this purpose, let us define two sequences of positive we define the generalized subsampling estimator of the σ ij , which is based on block lengths belonging to the set {b The idea behind the construction of the estimator ( 5) is as follows.Take any sequence b of integers such that b 1 ≤ b ≤ b 2 .Let t * be a random variable with uniform distribution on the set {1, 2, . . ., n − b + 1}.Note that the elements of the covariance matrix of the distribution τ b( θ n, b,t * − θ n ) (in the resampling world) are of the form (4), that is, cov n,b is a mean of the subsampling estimators (4) based on different block lengths.
n,b is also a natural candidate for the estimator.However, for the sake of simplicity, we decided to consider only σ (i,j) n,b in this article.
Let A = [a ij ] r×r be a known complex-valued matrix and let be a parameter of interests, which is a linear combination of elements of the matrix .Based on (5), we propose the following subsampling estimator of σ : Note that the proposed estimator is based on a standard subsampling estimator (4), which is averaged over different block lengths (one may also consider only single block length).We decided to provide very general formula for the estimator to make it suitable for variety of the possible future applications.Unfortunately, in contrary to the stationary case for any resampling method (including subsamping) for nonstationary processes (like considered in this article periodic and almost periodic time series) there is no result concerning the optimal subsample or block length choice.Many such results in the stationary setting are obtained using the mean square error expansions of the bias and the variance of the estimator.However, even in the simplest problem of the overall mean estimation for PC/APC time series, such expansion has not been obtained.There is also no heuristic approach for any standard parameter of interest.A broad simulation study was performed in Dudek and Potorski (2020) to check if known heuristic block/subsample length approaches for stationary time series may be adapted for periodic time series and none of them seemed to work.
In stationary case one may also find very interesting results (see, e.g., Politis, Romano, and Wolf 1999, chap. 10) stating that taking a linear combination of subsampling distribution estimators with different subsample sizes allows to construct an improved estimator in the sense of the second-order accuracy.
It is not possible to state, if such result may hold in our nonstationary case, since no Edgeworth expansion exists for PC/APC processes.Below we formulate conditions under which the subsampling estimator σ n is mean square consistent.
Theorem 3.1.Assume that {X t , t ∈ Z} is an α-mixing time series.Moreover, let b 1 and b 2 be subsample lengths such that 0 (i') OR there exists σ such that for any sequence of integers b (iii) for any 1 ≤ i ≤ r and any sequence b such that b 1 ≤ b ≤ b 2 we have convergence Then the estimator σ n (given by ( 7)) is mean square consistent, that is,

Spectral Density Estimation for General APC Time Series
In this section, we show how to apply results from the previous section in the problem of spectral density estimation of an APC time series with an almost periodic mean function, under general assumption that the set is unknown and finite.At first we formulate the assumptions and introduce some additional notation.Next, we provide an auxiliary central limit theorem for the Fourier coefficients of the mean function (Section 3.2.1).Additionally, we present the main results concerning consistency of our subsampling estimators (Sections 3.2.2 and 3.2.5).In Section 3.2.
(10) However, our main goal is to provide the exact form of the asymptotic covariance matrix.It turns out that its elements are linear combinations of the spectral density function.In a consecutive subsection, we use this feature to construct the estimator of P(ν, ω).
Our idea to obtain an estimate of P(ν, ω) is to apply the subsampling method in order to approximate the covariance matrix of the distribution N 4 (0, (ν, ω)).Then, having estimates of elements of the matrix (ν, ω) and using the formulas (11), we can easily get the estimate of P(ν, ω).In this way, in our approach we do not need to know the form of the set , that is, we do not need to know the Fourier frequencies of the mean function.
It is worth to indicate that in the literature dedicated to PC and APC processes, asymptotic normality results corresponding to Theorem 3.2, are treated as the preliminary step for construction of confidence intervals for parameters of interest.Since the asymptotic covariance matrix usually is of complicated form and depends on unknown parameters, in practice it is very difficult to estimate.Thus, in such situations to approximate quantiles of the asymptotic distribution, resampling methods are often applied (see, e.g., Dudek et al. 2014;Lenart 2016;Dudek 2018 and references therein).Although, in our problem we need to estimate the asymptotic covariance matrix (ν, ω) and not m(•, •), we also use a resampling method.More precisely, we adopt the idea of the subsampling estimator of the limiting variance examined in the stationary case in Carlstein (1986) and generalized to the nonstationary cases in Fukuchi (1999) and Politis, Romano, and Wolf (1999).It is necessary to indicate that the consistency of subsampling is not a sufficient condition to show that the subsampling estimator of the limiting variance (or covariance matrix) converges to the asymptotic variance (or the covariance matrix) and hence the mentioned consistency results are not sufficient for our purpose.Thus, in the next subsection, we introduce the consistent subsampling estimators of the linear combinations of the elements of the matrix (ν, ω) (see Theorem 3.2).To show their consistency, we adopt ideas presented in Fukuchi (1999), Carlstein (1986), and Politis, Romano, and Wolf (1999).
Theorem 3.3.Assume that the conditions of Theorem 3.2 hold.Additionally, assume that there exists ξ > 0 such that sup t X t 4+3ξ < ∞ and ∞ k=1 (k + 1) c−2 α ξ c+ξ (k) < ∞, where c is the smallest even integer such that c ≥ 4 + 2ξ .Then the conclusion of Theorem 3.1 holds for any complex valued matrix A = [a ij ] i,j=1...,4 and σ = σ (ν, ω) given by ( 6) with Setting a 13 = 1, a 24 = 1, a 14 = −i, a 23 = i and a ij = 0 otherwise, we get the estimate of 2π P(ν, ω) (see ( 11)).From now on the estimator σ n (ν, ω) obtained for this particular choice of coefficients a ij will be denoted by , that is, estimators calculated on the basis of one fixed subsample length b.Thus, we have where In the next section, we provide another interpretation of W sub n,b (ν, ω) and discuss its relationship with the classical spectral density estimator.

Another Look on Subsampling Spectral Density Estimator W sub n,b 1 ,b 2
In this section, we show that W sub n,b 1 ,b 2 (ν, ω) is in fact an average of generalized type of a periodogram-subsampling estimator for the spectral density function.Moreover, we present its relation to the classical spectral density estimator with the Bartlett window.
For the sake of clarity at first we recall a periodogramsubsampling estimator for the spectral density in a stationary case.
Let {X t , t ∈ Z} be a real-valued stationary time series with constant mean μ and spectral density function f (•) and (X 1 , . . ., X n ) be the available sample.Moreover, let where ) is the discrete Fourier transform of the sample.The subsampling version of (15) based on the subsample (Z t , . . ., Z t+b−1 ), t = 1, . . ., n − b + 1 is defined as follows: where d n,b,t (ω) = t+b−1 j=t Z j exp(−ijω), ω ∈ [0, 2π).Then, the subsampling periodogram estimator of the spectral density f is given by In ( 17), each subsample is demeaned using the estimator of the overall mean, however it is also possible to demean each subsample using the subsample average.In such a case at ω = 0, the estimator 2π f sub n,b (0) is the standard subsampling estimator of the variance of the sample mean or of the scaled spectral density at 0, 2π f (0) = lim n→∞ nvar(1/n n i=1 X i ) (see, e.g., Politis, Romano, and Wolf 1999;Politis and Romano 1994).Under appropriate mixing and moment assumptions, the estimator f sub n,b (ω) can be proven to estimate consistently f (ω).This will also hold true for some classes of nonstationary time series.In particular one may expect consistency for zero mean periodically and almost periodically correlated time series.The proof of this fact can be obtained, for example, following ideas of Soedjak (2002).Soedjak constructed an estimator using sequence of dependent sample paths of a zero mean strongly harmonizable process.Applying subsampling allows to replace the replicates of the process by the subsamples.
It seems natural that subsampling will work also for the time series with nonstationary mean structures (e.g., periodic, almost periodic) where there is a long-run pattern in the expected value of sample averages.Then instead of demeaning the original data, one may use subsampling to remove the mean at the level of the subsample averages.In a case of time series having mean structures with convergent long-run averages, this is the important feature for the validity of the method.And in fact this feature of subsampling is used by the estimator Note that in a nonzero mean case the estimator W sub n,b (ν, ω) introduced in this article, can be expressed as follows: where q = n − b + 1 and is the generalization to nonzero mean case of subsampling version of the dual-frequency periodogram.The dual-frequency periodogram was considered in zero mean case, for example, in Gorrostieta, Ombao, and von Sachs (2019) and Lenart (2011).According to our knowledge in the literature, there are no results for I n,b,t (ν, ω) when the set is unknown.
The subsampling estimator ( 18) is related with a classical estimator of the spectral density function.More specifically, as we see later W sub n,b (ν, ω) given by ( 14) (or equivalently by ( 18)) can be expressed using the classical spectral density estimator with the Bartlett window applied to the centered data Y t = X t − μ(t), that is, Lenart (2011) with L n = b, p. 293).Note that if the mean function μ(t) is unknown, then H n,b (ν, ω) cannot be used, while one can apply the estimator W sub n,b (ν, ω).Theorem below states relation between H n,b (ν, ω) and W sub n,b (ν, ω).
Theorem 3.4.Assume that conditions A1-A3 hold and that sup given by ( 14) can be expressed as follows: is deterministic and H n,b (ν, ω) is given by ( 19).
(a) (20) also holds under conditions of Theorem 3.2 with an additional assumption on the subsample length, that is, n = O(b 2+κ ) for some κ > 0, as n −→ ∞.(b) A relation analogous to (20) for the stationary case at zero frequency can be found in Politis, Romano, and Wolf (1999) (see Section 3.8.2).It was established for demeaned data, that is, under assumption that the set is known.Note also that since we assume that is unknown, H n,b (ν, ω) is not an estimator.
Note that for any fixed Also note that the set ˜ depends only of the set , that is, of the frequencies of the mean function.
The following corollary illustrates the relation between the limiting distributions of the rescaled W sub n,b (ν, ω) and H n,b (ν, ω).
Corollary 3.1.Assume that conditions of Theorem 3.4 hold.Denote and Then for any fixed (ν, ω) ∈ [0, 2π) 2 (i) if n/b 3 −→ 0 as n −→ ∞, then T 1n and T 2n have the same limiting distributions (assuming that T 2n is convergent in distribution); (ii) if n/b 3 → C 1 as n −→ ∞, where C 1 is some positive constant or infinity, then -for any (ν, ω) ∈ S, T 1n and T 2n have the same limiting distributions (assuming that T 2n is convergent in distribution); -forany(ν, ω) ∈ S the deterministic term n/b C n,b (ν, ω) is directly related to the mean of the distribution of T 1n .
Since n/b C n,b (ν, ω) does not have limit as n → ∞, the limiting distribution of T 1n does not exist.

Simulation Data Examples
Example below illustrates performance of our estimator Example 3.1.We consider a PC series X t = 5 sin (2π t/30) + 2 (sin (2π t/12) + 21/20) t , where t is a normalized Gaussian white noise.Thus, {X t } is composition of the PC white noise 2 sin 2π t 12 + 21 20 ε t and the periodic deterministic component 5 sin (2π t/30)), with ⊂ {λ k = 2kπ/12 : k = 0, . . ., 11}.It is easy to show that the spectral densities (denoted here by g λ k (•) = g k (•)) are constant on the interval [0, 2π) and g 0 = 6.41, g 1 = g 11 = −4.2i,g 2 = g 10 = −1, g 3 = • • • = g 9 = 0 (see Hurd and Miamee 2007, p. 170).To construct our estimator of spectral density function, we generated two samples of sizes n = 1000 and n=10,000, respectively.We set b , where x is the largest integer that is less than or equal to x.The simulation results are presented in Figures 1 and 2. Since the main conclusions do not depend of the chosen values of parameters, we decided to restrict the number of figures only to those that are necessary to illustrate performance of our estimator.
In Figure 1, one may observe the estimated rescaled spectral density 2π g 0 for different values of b.Let us recall that the main advantage of our approach is the ability to estimate the spectral densities without need of removing the mean function from the data.In this example the mean function is periodic and its true significant frequency ψ 0 is equal to 2π/30.First of all, one may notice that independently of the subsample length b our estimator provides an accurate estimate of the rescaled spectral density value at ψ 0 .However, in the neighborhood of ψ 0 , there are always two side-lobes representing the spectral leakage in the frequency domain.The size of lobes depends on the value of b, that is, on the subsample length.In Figure 2, we present the estimate of the spectral density functions on the bifrequency square.One may note that the real part of the estimate is concentrated on the main diagonal, while the imaginary part on two lines that are parallel to diagonal.Moreover, for the real and the imaginary parts, we observe leakage effect around vertical lines ν = ψ 0 , ν = 2π − ψ 0 and horizontal lines ω = ψ 0 , ω = 2π − ψ 0 .
As we mentioned above, the observed leakage is a consequence of applying the rectangular window of the length b to the data.In the result we do not calculate just the discrete Fourier transform (DFT) of {X t } but of the product of our signal {X t } and the rectangular window function.The Fourier transform of the product is equal to the convolution of the Fourier transforms.Moreover, the Fourier transform of the rectangular window function is the Dirichlet kernel function (see, e.g., Brillinger 2001, formula (3.3.5) and tab. 3.3.1 ).Finally, the convolution of the DFT of the signal with the Dirichlet kernel produces the observed side-lobes.
The leakage phenomenon is a typical problem in signal processing.Taking sample of the finite size for Fourier transform is equivalent to application of a rectangular window on a data stream.In the literature one may find many methods designed to reduce the leakage effect.The most popular are based on nonrectangular windows.However, reducing the spectral leakage results in decreasing the spectral resolution.More details on that topic can be found for instance in the book Proakis and Manolakis (2009).
Example 3.2.This example shows some advantages of our estimator over the classical estimators of the spectral density function.In contrary to Example 3.1 we considered time series with nonconstant spectral density functions.These are two AR(2) signals of the following form: where t is normalized Gaussian white noise and μ(t) = 1.6 sin(π t/3).Model M1 is the example of an AR(2) type stationary time series.Model M2 is having periodic mean function.Both models have the same spectral density function which concentrates the spectral mass around frequency π/3.We chose those two simple series to illustrate problems of classical estimator of spectral density function when the mean function is nonzero and to show that our method provides comparable results to the classical approach in the standard stationary series.
We estimated the spectral density function using our estimator W sub n,b 1 ,b 2 (ν, ω) (see Section 3.2.2) and the classical estimator based on smoothing the covariance function using the Bartlett window.We considered two sample sizes n = 500 and n = 1000 and block lengths b The length of smoothing window c was set c = (b 1 + b 2 )/2 .Since the conclusions do not depend of the sample size we decided to present results only for n = 1000 (see Figure 3).In panels on the left hand-side one may observe the generated realizations for models M1 and M2.Panels on the right-hand side contain the theoretical spectral density function (dashed line) together with its estimates obtained using our non-modified subsampling estimator (solid line) and the classical estimator with Bartlett window of length c (shadowed area).
In the case of M1, the mean function is equal to zero and hence assumption of the classical approach is fulfilled.As a result performances of the classical estimator and the subsampling estimator are comparable and obtained estimates are close to the true values.Model M2 is having periodic mean function with the true significant frequency equal to π/3.Note that the subsampling estimator provides estimate close to the true value of the spectral density function at the frequency π/3.At the same time the classical estimator fails at this frequency.Moreover, for M2 we observe the spectral leakage around π/3.
Examples above show that our estimator performs correctly and sometimes outperforms the classical one.In contrary to the classical approach, we do not need to remove the mean function from the data to be able to estimate the spectral density functions.That is very important in the real data applications since the period length of the mean function may not be known prior to the performed analysis on the data or the mean function may be almost periodic.In such cases, removing the mean function component from the data may be very challenging.
It should be emphasized that the leakage phenomenon that we observed appears naturally when DFT is applied to the data.We are aware that there exist variety of approaches trying to deal with the leakage problem.The aim of this article is not to provide the new solution for the leakage problem, however in the next section we propose the modified subsampling estimator that allows for significant reduction of the leakage effect.

Modified Subsampling Estimator
Below we introduce the modified subsampling estimator that significantly reduces the leakage problem.Instead of using , where v n is a deterministic sequence of real numbers such that v n → 0 as n → ∞.In other words, we consider a small neighborhood of (ν, ω).At first, we show mean square convergence of In the second step to reduce leakage effect, we use a finite number of estimators Let k > 0 be a fixed odd integer and ν, ω ∈ [0, 2π) be fixed frequencies.Moreover, let sequences v n,1 , . . ., v n,k be of the form v n,i = D i b −ρ i 1 , where 0 < ρ i < 1/4 and D i ∈ R, for i = 1, 2, . . ., k.To reduce the leakage effect we will use the  where the symbol "med" denotes the median.Applying M sub n,b 1 ,b 2 (ν, ω) to the points (ν, ω) that are in the area of the side-lobes, especially in the case when the side-lobes are steep, allows to reduce size of the lobes.On the other hand, when (ν, ω) is in the area, in which values of the variance estimator are not changing much (i.e., it is not in the area of the side-lobes), then using M sub In Figure 4, one may observe the estimated rescaled spectral density 2π g 0 for different values of b.Definitely, using the estimator M sub n,b 1 ,b 2 allowed for significant reduction of the leakage effect.The side-lobes are no longer present around λ 0 = 2π/30.The same conclusion can be drawn looking on the estimates on the bifrequency square (see Figure 5).

Conclusions
In this article, we proposed a new nonparametric approach for spectral density estimation, which does not require demeaning of the data and is based on the subsampling method.Our estimator was constructed for almost periodically correlated sequences, which are an important subclass of nonstationary time series.Additionally, we introduced the modified estimator, which allows to reduce the leakage effect.We showed the ability of our approach to analyze the cyclical fluctuations and to determine the period length on the examples of typical economic datasets.Our method may be used to construct testing tools to recognize the type of cyclical fluctuations, that is, to decide if they are described by the periodic mean function or by the higher-order moments (for more details and the economic data examples we refer the reader to the supplementary materials).
Presented results are providing new open problems concerning subsampling estimation for the higher-order moments.Since the classical approach for the higher-order spectral estimation also requires removing the mean function from the data, we believe that ideas presented in this article can be adapted in higher-order spectral estimation problem under nonzero mean assumption.This and other related issues are the subject of the current research of the authors.
Additionally, assume that (i) there exists σ such that for any sequence of integers b = b(n) such that b 1 ≤ b ≤ b 2 , and any sequence of integers
ω) and at the same time problem of leakage, which inspired us to propose in the next subsection a modification of W sub n,b 1 ,b 2 (ν, ω).

Figure 3 .
Figure 3. From top results for models M1-M2, respectively.Left panels: generated realizations; right panels: classical spectral density estimates obtained using Bartlett window of length c (shadowed area), subsampling estimates (solid line), theoretical spectral density functions (dashed line).
define the new estimator as their median.Theorem 3.5.Assume that the conditions of Theorem 3.3 hold.Additionally, assume that v n = Db −ρ 1 , where 0 < ρ < 1/4 and D ∈ R.Then, as n −→ ∞ we have (a) W sub n,b 1 ,b
3, we discuss the relation of our estimator with the periodogram-subsampling estimator of the spectral density function and with classical spectral density estimator based on the Bartlett window.Section 3.2.4provides short illustrative examples with simulated data.Let X c n , . . ., X c n +d n −1 be a sample of the size d n from an APC time series {X t , t ∈ Z} with the mean function μ(t), where {c n } n∈Z and {d n } n∈Z are arbitrary sequences of positive integers such that d n → ∞ as n → ∞.Let us recall that when the set is finite, the mean function μ(t) can be expressed as follows: μ(t) = t ∈ Z} is a real-valued APC time series such that the sets and are finite and unknown; A2 {X t , t ∈ Z} is α-mixing (for definition see, e.g.
t = μ(t) + Y t , where Y t is a zero mean harmonizable sequence (in the sense of Loève 1963; see Hurd and Miamee