Estimation of a nonparametric regression spectrum for multivariate time series

Statistical Methodology 6 (2009) 202–222

Contents lists available at ScienceDirect

Statistical Methodology journal homepage: www.elsevier.com/locate/stamet

Estimation of a nonparametric regression spectrum for multivariate time series Jan Beran a,∗ , Mark A. Heiler a,b a

Department of Mathematics and Statistics, University of Konstanz, Germany

b

Credit Suisse, Switzerland

article

info

Article history: Received 8 December 2007 Received in revised form 15 July 2008 Accepted 6 September 2008 Keywords: Periodogram Cross-spectrum Regression spectrum Phase Wavelets

a b s t r a c t Estimation of a nonparametric regression spectrum based on the periodogram is considered. Neither trend estimation nor smoothing of the periodogram is required. Alternatively, for cases where spectral estimation of phase shifts fails and the shift does not depend on frequency, a time domain estimator of the lag-shift is defined. Asymptotic properties of the frequency and time domain estimators are derived. Simulations and a data example illustrate the methods. © 2008 Elsevier B.V. All rights reserved.

1. Introduction Consider a multivariate time series Y(i) = (Y1 (i), . . . , Yp (i))T of the form Y(i) = f(ti ) + (i),

(1)

where ti = i/n (i = 1, . . . , n), f(t ) = (f1 (t ), . . . , fp (t )) ∈ C (t ∈ R) is a multivariate deterministic trend function and (i) = (ε1 (i), . . . , εp (i))T is a zero mean stationary process. In this model, dependence between two components Yr and Ys can occur due to two reasons: 1. dependence between r and s , and 2. dependence due to similarities in the underlying deterministic components fr and fs . In the first case, linear dependence is characterized by cross-correlations, the cross-spectrum, coherency and phase shift between r and s (see e.g. standard books such as [1,2]). For the second case, Beran and Heiler [3] introduced a nonparametric regression cross-spectrum. In the present paper, we consider estimation of the regression cross-spectrum based on the periodogram, and frequency T

∗

Corresponding author. E-mail address: [email protected] (J. Beran).

1572-3127/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.stamet.2008.09.001

p

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

203

and time domain estimation of possible phase shifts. Fig. 5 shows a typical example where the nonparametric regression spectrum leads to interesting insights. The bivariate series consists of the Southern Oscillation Index (Fig. 5a) and recruitment of new fish in the central Pacific Ocean (Fig. 5b), ranging from 1950 to 1987 over a period of n = 453 months [4]. Both series have strong deterministic components that are related to each other. The analysis in Section 5.2 shows that there are two levels of dependencies, namely between the long-term trends (El Niño effect) of both series, and between the deterministic seasonal components respectively. The paper is organized as follows. Basic definitions from [3] are summarized briefly in Section 2. An estimator of the regression spectrum, its modulus and the phase spectrum, based on the periodogram, is discussed in Section 3, together with asymptotic properties. The asymptotic results imply in particular that phase estimates can be highly unreliable for frequencies with low amplitude spectrum. In fact, examples in Section 5 illustrate that estimation of time delays from the raw plot of the (estimated) regression phase spectrum is virtually impossible. The problem is resolved by applying an algorithm that downweighs or eliminates unreliable frequencies. For cases where the number of relevant frequencies is too small, an alternative procedure for estimating time shifts between trend functions is presented in Section 4. Simulations and a data example in Section 5 illustrate the methods. Proofs are given in the Appendix. 2. Definition of the regression cross-covariance and spectrum Under suitable regularity assumptions on f, [3] define the regression (cross-)covariance function

0(u) = [γrs (u)]r ,s=1,...,p and the regression (cross-)correlation function R(u) = [ρrs (u)]r ,s=1,...,p of f(t ) by

γrs (u) = hfr (· + u), fs (·)i =

1

Z

fr (t + u)fs (t )dt 0

and

γrs (u) (2) ρrs (u) = √ γrr (0)γss (0) R1 where o f (x)dx is assumed to be zero and f(1 + u) = f(u) (0 < u ≤ 1). Note that here, trend components that cannot be extended periodically beyond t = 1 are assumed to have been removed, or to be negligible. For t ∈ [0, 1] and fr ∈ L2 [0, 1], we have ∞ X f(t ) = a(j)ei2π jt , j=−∞

where a(j) = (a1 (j), a2 (j), . . . , ap (j))T ∈ Cp are given by ar (j) = hfr , ei2π j· i =

1

Z

fr (t )e−i2π jt dt .

0

Hence, the regression spectrum at frequency j is defined as the sequence of p × p-matrices H(j) = [hrs (j)]r ,s=1,...,p (j ∈ Z) with H(j) = a(j)aT (j). The regression spectrum and covariance function are closely linked by

0(u) =

∞ X

H(j)ei2π ju

j=−∞

and H(j) =

Z

1 2

− 21

e−i2π ju 0(u)du.

204

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

Using polar representation of H(j),

|ar (j)as (j)| hrs (j) = √ exp(iφrs (j)) γrr (0) · γss (0) γrr (0) · γss (0)

h˜ rs (j) = √

is called the standardized regression spectrum of f,

|hrs (j)| |hrs (j)| κrs (j) = √ = kf r k · k f s k γrr (0) · γss (0) = rP

|ar (j)as (j)| P |ar (l)|2 |as (m)|2 m

l

is the relative spectral modulus and φrs (j) the phase shift at frequency j. The relative spectral modulus (or coherence) κrs (j) defined above can assume any number between 0 and 1, thus giving a relative measure of the contribution of frequency j to the cross-covariance. Remark 1. If two components fr and fs are shifted versions of each other, i.e. fs (t ) = cfr (t + ∆) for some ∆, c 6= 0, then

|ar (j)|2 κrs (j) = P |ar (l)|2

(3)

l

and the phase shift

φrs (j) = −2π j∆ is a linear function of the shift parameter ∆. 3. Estimation in the spectral domain 3.1. The periodogram In practice, the trend component f is usually unknown. Beran and Heiler [3] propose an estimator of the cross-spectrum based on a trend estimator ˆf obtained by wavelet thresholding. In this section, we consider a direct estimator of the regression spectrum based on the periodogram. This has two main advantages. First of all, the estimator is simple and does not require nonparametric estimation of the trend function f. The second advantage is that estimated values at different frequencies are asymptotically independent. Given n observations of a multivariate vector Y(i) (i = 1, . . . , n), the periodogram of Y(i) at frequency ωj = 2π j/n, ωj ∈ [−π , π], is defined by I(ωj ) =

1

n X

n

! Y(s) exp(−iωj s)

s=1

Moreover, let A(ωj ) =

X

f(tk ) exp(−iωj k)

B(ωj ) =

X

(k) exp(−iωj k),

and

n X t =1

!T Y(t ) exp(iωj t )

.

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

205

respectively. Then

  1 T I (ωj ) = I;rs (ωj ) 1≤r ,s≤p = B(ωj )B(ωj ) , n is the periodogram of the multivariate stationary series (i) = (1 (i), . . . , p (i))T . The deterministic counterpart,

  1 T If (ωj ) = If ;rs (ωj ) 1≤r ,s≤p = A(ωj )A(ωj ) n will be called regression periodogram of f. It can be seen by straightforward calculations that under model (1) with non-constant f the diagonal elements of the periodogram I(ωj ) are of the order O(n) and are dominated by If (ωj ). This is the essential reason why the regression spectrum can be estimated directly from the periodogram. Specific results on the asymptotic distribution are given in the following two theorems. Theorem 1. Denote by H = [hrs ]r ,s=1,...,p the regression spectrum of f, and suppose that (i) are independent, identically distributed zero mean random vectors with non-singular covariance matrix 6 = (Σrs )1≤r ,s≤p , and existing fourth moments. Then, for each pair (r , s), and Fourier frequencies 0 ≤ ωj1 < · · · < ωjk ≤ π , (k ∈ N), the following hold. (i) n−1 Irs (ωj ) − hrs (j) = Op (n−1/2 ),   (ii) E n−1 Irs (ωj ) − hrs (j) = O(n−1 ), (iii)

√

n[n−1 Irs (ωj1 ) − hrs (j1 ), . . . , n−1 Irs (ωjk ) − hrs (jk )]T converges in distribution to a k-dimensional normal random vector with mean 0 and covariances lim n cov(n−1 Irs (ωj ), n−1 Irs (ωj0 )) = 0,

n→∞

for ωj 6= ωj0 , and lim n var(n−1 Irs (ωj )) = Σrr |as (j)|2 + Σss |ar (j)|2 .

n→∞

(4)

Remark 2. More specifically, the proof of the theorem implies n cov(n−1 Irs (ωj ), n−1 Irs (ωj0 )) = O(n−2 ), for ωj 6= ωj0 , and n var(n−1 Irs (ωj )) = Σrr |as (j)|2 + Σss |ar (j)|2 + Rn

(5)

with Σrs denoting the (r , s)th entry of the matrix 6 and Rn = n−1 Σrr Σss + O(n−2 ) for 0 < ωj < π , and Rn = n−1 (Σrr Σss + Σrs Σsr ) + O(n−2 ) for ωj ∈ {0, π}. Remark 3. The periodogram I(ω) contains a stochastic and a deterministic part. This carries over to the asymptotic variance at frequency j. The main component is given by a product of the variance of the stationary part and the regression spectrum. The variance in Theorem 1 is of order O(n−1 ) so that n−1 Irs (ωj ) is an asymptotically consistent estimator of hrs (j). This is in contrast to spectral estimation for stationary processes where the periodogram needs to be smoothed. Note also that, in contrast to the wavelet estimator in [3], the estimators ˆf(ω) = n−1 I(ω) at different frequencies are asymptotically independent. This facilitates estimation of the coherence and phase spectrum (see results below).

206

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

Theorem 1 can now be extended to linear error processes. Thus, assume

(i) =

∞ X

A(j)Z(i − j),

j=−∞

where A(j) = (Alk (j))1≤l,k≤p are p × p-matrices such that for all pairs (l, k),

X

1

|Alk (j)||j| 2 < ∞,

j∈Z

and Z(i) are independent identically distributed with zero mean and non-singular covariance matrix 6. Also, denote by





h = h;r ,s r ,s=1,2,...,p the cross-spectral density of (i). Theorem 1 can be generalized to: Theorem 2. Let (i) be a linear process as defined above. Then, with the same notation as in Theorem 1, (i) n−1 Irs (ωj ) − hrs (j) = Op (n−1/2 ); (ii) E n−1 Irs (ωj ) − hrs (j) = O(n−1 ); √ (iii) n[n−1 Irs (ωj1 ) − hrs (j1 ), . . . , n−1 Irs (ωjk ) − hrs (jk )]T converges in distribution to a k-dimensional normal random variable with mean zero, lim n cov(n−1 Irs (ωj ), n−1 Irs (ωj0 )) = 0

n→∞

for ωj 6= ωj0 , and lim n var(n−1 Irs (ωj )) = 2π h;rr (ωj )|as (j)|2 + 2π h;ss (ωj )|ar (j)|2 .

n→∞

Remark 4. More specifically, the proof of the theorem implies n cov(n−1 Irs (ωj ), n−1 Irs (ωj0 )) = O(n−2 ) for ωj 6= ωj0 , and n var(n−1 Irs (ωj )) = 2π h;rr (ωj )|as (j)|2 + 2π h;ss (ωj )|ar (j)|2 + Rn

(6)

with Rn = n−1 (2π )2 h;rr (ωj )h;ss (ωj ) + O(n−2 ) for 0 < ωj < π , and

  Rn = n−1 (2π )2 h;rr (ωj )h;ss (ωj ) + h;rs (ωj )h;sr (ωj ) for ωj ∈ {0, π}. Remark 5. Recalling that hrs (j) = ar (j)as (j), we see that the asymptotic variance of Irs (ωj ) is equal to 2π times the product of the regression and the stationary spectrum. Remark 6. The result in Theorem 2 can be generalized to pairs (r , s) and (r 0 , s0 ), 1 ≤ r , s, r 0 , s0 ≤ p. For Fourier frequencies ωj , ωj0 we have lim n cov(n−1 Irs (ωj ), n−1 Ir 0 s0 (ωj0 )) = O(n−2 )

n→∞

for j 6= j0 and lim n cov(n−1 Irs (ωj ), n−1 Ir 0 s0 (ωj0 )) = 2π hss0 (j)h;rr 0 (ωj ) + 2π hrr 0 (j)h;ss0 (ωj ) + O(n−1 )

n→∞

for j = j0 .

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

207

3.2. Estimation of modulus, phase spectrum and phase shift We now consider estimation of the modulus and the phase spectrum. As a special case, the asymptotic results will be applied to the estimation of a constant phase shift. Denote by crs (j) = Re{hrs (j)} and qrs (j) = −Im{hrs (j)} the real part and the imaginary part with reversed sign of hrs (j) = ar (j)as (j) respectively. Thus, hrs (j) = crs (j) − iqrs (j). Similarly, we define for the cross-spectral density of ε(i) the quantities cε;rs (j) = Re{hε;rs (j)} and qε;rs (j) = −Im{hε;rs (j)}. The joint asymptotic distribution of cˆrs (j), qˆ rs (j) follows directly from Theorem 1 and Remark 6. Lemma 1. Define the estimators 1 h

cˆrs (j) =

2n

Irs (ωj ) + Irs (ωj )

i

and 1 h

qˆ rs (j) = −

2ni

i

Irs (ωj ) − Irs (ωj ) .

Then,

ζn =

√ 

n cˆrs (j) − crs (j), qˆ rs (j) − qrs (j)

T

converges in distribution   to a bivariate normal zero mean random variable with asymptotic covariance matrix M (j) = Mij (j) i,j=1,2 given by M11 (j) = lim n var(ˆcrs (j))

n→∞     = π hss (j)h;rr (ωj ) + hrr (j)h;ss (ωj ) + 2π crs (j)c;rs (ωj ) − qrs (j)q;rs (ωj ) ,

M22 (j) = lim n var(ˆqrs (j))

n→∞     = π hss (j)h;rr (ωj ) + hrr (j)h;ss (ωj ) − 2π crs (j)c;rs (ωj ) − qrs (j)q;rs (ωj ) ,

and M12 (j) = lim n cov(ˆcrs (j), qˆ rs (j)) n→∞

  = 2π qrs (j)c;rs (ωj ) + crs (j)q;rs (ωj ) . Based on Lemma 1 we may now construct consistent estimators of the spectral modulus and the phase shift. Using the notation hrs (j) = κrs∗ (j) exp(iφrs (j)) estimators of the spectral modulus and phase shift respectively are defined by

κˆ rs∗ (j) =

q

cˆrs2 (j) + qˆ 2rs (j)

(7)

208

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

and by

  qˆ rs (j) φˆ rs (j) = arg Irs (ωj ) = arctan − . cˆrs (j) The asymptotic distribution of κˆ rs∗ (j) and φˆ rs (j) is given in the following corollaries.

(8)

Corollary 1. Let κrs∗ (j) > 0 and let the assumptions of Theorem 2 hold, then

κˆ rs∗ (j) − κrs∗ (j) = Op (n−1/2 ) uniformly for all j ∈ Z. Furthermore,

√ 

d

2 n κˆ rs∗ (j) − κrs∗ (j) → N (0, σκ; rs (j)),



where crs2 (j)M11 (j) + q2rs (j)M22 (j) + 2crs (j)qrs (j)M12 (j)

2 σκ; rs (j) =

 2 κrs∗ (j)

.

(9)

Moreover, lim n cov(κˆ rs∗ (j), κˆ rs∗ (j0 )) = O(n−2 )

n→∞

(j 6= j0 ).

Corollary 2. Let κrs∗ (j) > 0 with φrs (j) = arg hrs (j) and φˆ rs (j) as above. Then,

φˆ rs (j) − φrs (j) = Op (n−1/2 ) uniformly for all j ∈ Z. Furthermore,

√ h

i

d

2 n φˆ rs (j) − φrs (j) → N (0, σφ; rs (j)),

where q2rs (j)M11 (j) + crs2 (j)M22 (j) − 2crs (j)qrs (j)M12 (j)

2 σφ; rs (j) =

 4 κrs∗ (j)

.

Moreover, lim n cov(φˆ rs (j), φˆ rs (j0 )) = O(n−2 )

n→∞

(j 6= j0 ).

Note that the variance of the phase spectrum is small whenever the modulus κrs∗ (j) is large and vice versa. Accurate estimation of the phase spectrum may therefore only be expected for frequencies j where the amplitude spectrum is large. Examples in Section 5 illustrate that often most frequencies have to be omitted in the estimation of phase shifts. For instance, in the case of a simple shift between fr and fs , i.e. fs (t ) = cfr (t + ∆) and hence

φrs (j) = −2π j∆, the following algorithm can be applied: 1. Calculate ˆf(ωj ) = n−1 I(ωj ); 2. Define J∗ =



j : κˆ rs∗ (j) > c ·

q

2 σκ; rs (j)



for a suitably chosen c ∈ R; 3. Estimate the phase shift by applying a local robust regression to the points {(j, φrs (j)) : j ∈ J ∗ }, taking into account possible jumps modulo 2π .

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

209

4. Lag estimation in the time domain In the previous sections we considered estimation of the regression cross-spectrum based on the periodogram and derived a method for estimating lead–lag effects in the trend components. The proposed algorithm is based on a set of significant common frequencies that can be used for estimating the slope of the phase line. Problems with this algorithm are expected if the set of common frequencies is too small to identify a slope in the phase plot. This is the case, for instance, if the deterministic components have a Fourier series representation with a small number of harmonic components. Phase shifts may then have to be identified by examining regression cross-correlations instead of the phase spectrum. In the case of a simple shift that does not depend on frequency, the time delay between two trend components can be estimated by identifying the maximum of the cross-correlation. Thus, for each pair (r , s), 1 ≤ r , s ≤ p, denote the set of local maxima by

M = {u ∈ [−1, 1] : γrs0 (u) = 0, γrs00 (u) < 0} and umax = argmax{γrs (u) : u ∈ M} rs where γrs is the cross-auto-correlation defined in Section 1. An estimator of umax is then defined by rs uˆ max = argmax{γˆrs (u) : u ∈ Mˆ } rs where

Mˆ = {u ∈ [−1, 1] : γˆrs0 (u) = 0, γˆrs00 (u) < 0} and γˆrs is a suitable consistent estimator of γrs (u). More specifically, here γˆrs will be defined using a wavelet estimator of the trend function f ∈ Rp in the definition of γrs . Thus, given a father and a mother wavelet φ(·), ψ(·) ∈ L2 (R) and the corresponding wavelet basis l

φl,k (x) = 2 2 φ(2l x − k) and j

ψj,k (x) = 2 2 ψ(2j x − k),

k, j ∈ Z,

we define

γˆrs (u) =

1

Z

fˆr (t + u)fˆs (t )dt 0

with fˆr (t ) :=

X

αˆ l(,rk) φl,k (t ) +

k

Jn X X j ≥l

w ˆ j(,rk) βˆ j(,rk) ψj,k (t ),

(10)

k

where

αˆ l(,rk) =

n 1X

n u =1

φl,k (tu )Yr (u),

(11)

ψj,k (tu )Y (u),

(12)

and

βˆ j(,rk) =

n 1X

n u =1

210

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

(r )

(r )

for some Jn → ∞, and w ˆ j,k := 1{|βˆ j,k | ≥ (r )

q

(r )

var(βˆ j,k )λj }. For the choice of the threshold λj see e.g. [5– (r )

(r )

(r )

7], among others. Denote by α(r ) = {αl,k : l, k ∈ Z, αl,k 6= 0} and β(r ) = {βj,k : j ≥ l, k ∈ Z, βj,k 6= 0} (r = 1, . . . , p) the coefficients in the wavelet representation of the components of f, i.e. fr (t ) :=

X

αl(,rk) φl,k (t ) +

Jn X X j ≥l

k

βj(,rk) ψj,k (t ).

(13)

k

(r )

(r )

As in [5], we assume that for each r, the number of non-zero coefficients αl,k and βj,k is finite. Let r , s ∈ {1, 2, . . . , p} be fixed. Defining

θ0 = (α(r ) , β(r ) , α(s) , β(s) ) we may then write

γrs (u) = γrs (u, θ0 ) where γrs depends continuously on θo ,

M = M(θ0 ) = {u ∈ [−1, 1] : γrs0 (u, θ0 ) = 0, γrs00 (u, θ0 ) < 0} and umax = argmax{|γrs (u, θ0 )| : u ∈ M(θ0 )}. rs

(14)

The estimator of umax is then defined by rs

n



o

uˆ max = argmax γrs (u, θˆ ) : u ∈ M(θˆ ) . rs

(15)

To ensure existence, uniqueness and consistency of the estimator the following assumptions will be used. (A1) φ, ψ have compact support and are of finite variation; (A2) fr (r = 1, 2, . . . , p) are as defined in Section 1 and of bounded variation; (A3) The cumulants cm;r (u1 , . . . , um−1 ) = cum {εr (i + u1 ), . . . , εr (i + um−1 ), εr (i)} of εr (i) exist, are absolutely summable, i.e. Cm;r =

X

cm;r (u1 , . . . , um−1 ) < ∞.

u1 ,...,um−1

Moreover, εr has covariances γεr (k) such that ∞ X kγε (k) < ∞; r k=−∞

(A4) dim(θo ) < ∞; (A5) For z in a small neighborhood of 0,

X

Cm z m < ∞;

m

(A6) As n → ∞, we have Jn → ∞, n2−Jn /2 → ∞, 2j/2 λj = o(n1/2 ) (j = l, l + 1, . . . , Jn ) and Jn X

2j/2 exp(−2λ2j /(1 + η)) = o(1)

j >l

for some η > 0; (A7) γrs (u, θ ) is twice continuously differentiable with respect to u and θ ; (A8) |γrs (u, θ0 )| has a unique maximum at umax rs .

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

211

Asymptotic properties of uˆ max are given in the following theorem. rs Theorem 3. Under assumptions (A1)–(A8) we have, for 1 ≤ r , s ≤ p, uˆ max − umax = Op n−1/2 rs rs




and

√

d

n uˆ max − umax −→ N (0, τu,rs ), rs rs




with

τu,rs (θ0 ) =

1 2 (γrs00 (umax rs , θ0 ))



∂ 0 max γ (u , θ0 ) ∂θ rs rs

T

var(θˆ )



 ∂ 0 max γrs (urs , θ0 ) . ∂θ

Remark. A similar result for univariate stationary time series can be found in [12]. 5. Examples 5.1. Simulations Consider model (1) with f1 (x) a piecewise constant function as displayed in Fig. 1c, f2 (x) = f1 (x + ∆) with ∆ = .0625, 1 (i), 2 (i) independent and identically distributed N (0, σ2 ), σ2 = 9 and corr(1 (i), 2 (j)) = 0. A simulated sample path of Y(i) = (Y1 , Y2 )T (i = 1, 2, . . . , 2048) with Yj (i) = fj (ti ) + j (i) (j = 1, 2) is displayed in Fig. 1a and b. The regression amplitude and phase spectrum for these trend components are shown in Fig. 1e and f. Estimates of the regression amplitude and phase spectrum obtained from n−1 I are shown in Fig. 1g and h respectively. Fig. 1g, h illustrate that the common frequencies can be identified quite accurately in the amplitude spectrum, whereas the phase spectrum is heavily disturbed by the random noise components 1 , 2 . This is expected in view of Corollaries 1 and 2. It is therefore essential to use important common frequencies only, when estimating the regression phase spectrum. Fig. 2a through d display results of a small simulation study where the simulated and true variance of the amplitude spectrum according to (9) are compared for different sample sizes. At each frequency N = 200 simulations are carried√out and the amplitude spectrum was estimated. The empirical standard deviations multiplied by n are plotted in Fig. 2a through d (black line) together with their asymptotic counterparts (red line). Convergence to the asymptotic standard deviation is apparently faster for frequencies where the amplitude spectrum is large. These are exactly the ones which are used for estimating time shifts. Fig. 3a shows the amplitude spectrum and the phase estimate for one simulated series. Frequencies where the estimated amplitude spectrum is above four times its standard deviation are highlighted by black squares in Fig. 3a and b. The resulting estimated phase line in Fig. 3b is obtained by linear regression using these points only, taking into account jumps modulo 2π . The red lines indicate 99% confidence intervals for the regression slope. The regression line in Fig. 3b, with slope around 0.39, is obviously very similar to the true phase spectrum (Fig. 1f) with slope 0.30. For a more systematic ˆ , a small simulation study was carried out, with sample illustration of finite sample properties of ∆ ˆ (Fig. 4) sizes n = 512, 1024, 2048 and n = 4096. In each case we ran 200 simulations. Boxplots of ∆ based on 200 simulations (for each n) illustrate that estimation of ∆ is very difficult for n = 512. ˆ improves fast, however, with increasing sample size. A detailed summary of this The accuracy of ∆ simulation study is given in Table 1. Finally, we examine in how far confidence intervals for ∆, based on weighted linear regression ˆ (with weights and residual variances obtained from Corollary 2) have the desired coverage of (j, φ) probability. For each simulated series, the six frequencies with largest cross-spectral modulus were used in the regression, and 95% confidence intervals were calculated using estimated variances of φˆ at

212

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

Fig. 1. Simulated series of length n = 2048 (a, b), as defined in Section 5.1, trend components (c, d), and amplitude and phase of the true regression spectrum (e, f). Estimates of the amplitude and phase spectrum based on the periodogram are given in (g) and (h).

Table 1 Summary statistics of lag estimates. For each sample size, 200 simulations were carried out n

512

1024

2048

4096

True value Median Mean Std. dev.

0.0625 0.06012 0.05856 0.03816

0.0625 0.06073 0.05838 0.02275

0.0625 0.06322 0.06306 0.00918

0.0625 0.06300 0.06262 0.005528

these frequencies. The coverage percentages, based on 1000 simulations, turned out to be close to the desired values, namely 93.9%, 93.5%, 94.8% and 94.8% for n = 512, 1024, 2048 and 4096 respectively. Compared to the method in [3], the approach proposed here has the advantage of being much simpler, because no trend estimation is required. The question arises, however, whether using the raw data may not lead to less accurate estimates. A theoretical answer is not easy to obtain and, in general, depends on the signal-to-noise ratio and the type of trend function. We will not pursue this here in detail. A small simulation study illustrates, however, that the periodogram based

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

213

Fig. 2. Simulated and theoretical (red line) standard deviations for the amplitude spectrum for various sample sizes n = 512, 1024, 2048 and n = 4096. For each sample size, 200 simulations were carried out at each individual frequency. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Estimates of the amplitude and phase spectrum based on the periodogram for frequencies j = 1, . . . , 50. Important common frequencies are marked by black squares. The estimated phase line (solid line) is based on these frequencies. The red lines represent 99% confidence intervals. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

estimate can compete with the wavelet approach, at least for large enough samples. For n = 1024, 2048, 4096 and 8192, the following models were simulated: (a) Yj (i) = fj (ti ) + j (i) (j = 1, 2); (b) Yj (i) = fj (ti ) + 3j (i) (j = 1, 2) and (c) Yj (i) = gj (ti ) + j (i) (j = 1, 2) where f2 (t ) = f2 (t + ∆) with f1 and ∆ as given above and g2 (t ) = g1 (t + ∆) with



g1 (t ) = 1 t <

1 2



(t − 1) + 1



1 2

≤t≤

3 4





t +1 t >

3 4



(t − 2).

The comparison between models (a) and (b) illustrates the effect of an increased variance in the noise component, whereas models (a) and (c) differ with respect to the trend function. Simulated series of length n = 1024 are displayed in Fig. 5a and b. Fig. 6a through c show, for each n, boxplots of estimated values of ∆ based on the periodogram (left) and the wavelet (right) approach respectively. Note that here, a different sample size dependent threshold for including frequencies was used. It

214

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

Fig. 4. Boxplots of simulated lag estimates for n = 512, 1024, 2048 and n = 4096. For each sample size, 200 simulations were carried out. The true value 0.0625 is marked by a horizontal line.

Fig. 5. Simulated series of length n = 1024 for Models b (a) and c (b) respectively.

was chosen such that, empirically, it is close to optimal for the wavelet method. Nevertheless, the wavelet method is not uniformly better. For both methods, increasing the variance of the residual ˆ (Fig. 6a and b). With respect to the bias, the periodogram process leads to higher variability in ∆ based estimate is affected more whereas this is not the case for the variance. The trend function in ˆ . The reason may be that g1 is mostly linear so that model (c) also leads to increased variability of ∆ a shift is more difficult to detect than for a function with a large number of abrupt changes (such as f1 ). A difference between the periodogram and the wavelet approach is mainly visible for smaller P ˆ i − ∆i )2 versus n, in sample sizes. This is confirmed by plots of mean squared errors MSE = N −1 (∆ log–log-coordinates (Fig. 6d through f). For n = 8192, both methods essentially lead to the same mean squared errors. Only for model (b), the wavelet approach appears to be considerably better for smaller sample sizes. There, noise appears to be too large to rely on the raw data only. This effect disappears for n = 8192.

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

215

Fig. 6. Boxplots of simulated lag estimates (a–c) for n = 1024, 2048, 4096 and 8192. For each sample size, the left boxplot corresponds to the periodogram based estimate and the right boxplot to the wavelet based estimate (as proposed in [3]) respectively. The true value of ∆ = 0.0625 is marked by a horizontal line. (d–f) show simulated mean squared errors plotted against n in log–log-coordinates, for the periodogram (circles) and wavelet based (triangles) estimates and the three models respectively.

5.2. El Niño and recruitment of new fish Fig. 7a and b display the components of the bivariate time series consisting of the Southern Oscillation Index (SOI) and recruitment (amount) of new fish in the central Pacific Ocean (Fig. 7a and b), ranging from 1950 to 1987 over a period of n = 453 months. The SOI relates changes in air pressure to the temperature of the ocean at the surface. The data set can be found in [4]. Both time series exhibit cyclic components. The dominating periodic component in the SOI has a period of 12 months. The second series oscillates with a lower frequency, but a 12-month cycle is visible as well, at least in parts of the series. This is most obvious when looking at the amplitude of the estimated regression cross-spectrum in Fig. 7c which shows a dominating frequency at j = 38 indicating a period of 453/38 ≈ 12 months. In addition, a certain number of moderate contributions are present at low frequencies. The slight influence of low frequency components is also visible directly in the SOI series, in that the mean and variability seem to be changing slowly. This feature is often referred to in the literature as the El Niño effect. We now proceed as in the simulated example. The estimated phase line in Fig. 7f is based on frequencies where the amplitude spectrum is large. The corresponding points are marked by black squares in Fig. 7e, f. Focussing on these points only, one can detect a linear structure. We may thus assume that there is a frequency-independent shift ∆. Linear regression yields a slope ˆ = 0.1/(2π ) ≈ 0.016 which corresponds to 453 · 0.1/(2π ) ≈ 7.2 estimate of about 0.1, and thus ∆ months. This indicates that the SOI signal leads the recruitment of new fish by about seven to eight months. However, in view of Fig. 7c, better insight may be gained by separating high and low frequency components in the second series. We therefore carry out the analysis at separate levels of resolution. Fig. 8a and b show trend estimates fˆ1 (SOI) and fˆ2 (fish recruitment) obtained by wavelet thresholding with s20-wavelets. The second trend function is decomposed further by separating the three coarsest resolution levels (Fig. 8c) from the fourth level (Fig. 8d). The fourth level represents the 12-month cycle, whereas the coarser parts (levels one to three, D4-S6) represent four to five year cycles that may

216

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

Fig. 7. Monthly SOI and numbers of new fish (a, b) for the period 1950–1987, estimated amplitude of cross-spectrum (c, d), cross-amplitude spectrum (e) with line representing three times the standard deviation and phase spectrum (f). Important frequencies are marked in both plots by black squares. Estimation of the phase line is based on these frequencies. The red lines in (f) represent 99% confidence intervals for the regression line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

be associated with corresponding cycles in the warming of the Pacific Ocean. The estimated amplitude cross-spectrum with significant frequencies marked by black squares are displayed in Fig. 8e and the resulting phase line with corresponding confidence intervals is presented in Fig. 8f. One notices again the distinct linear structure over this particular set of frequencies. The slope of this line is given by 0.17 indicating a lead of SOI of 453 · 0.17/(2π ) ≈ 12 months. Since annual seasonality is the only common periodic oscillation between SOI and the component D3 (Fig. 8g), time delays between these two series are estimated in the time domain. The regression correlations between the SOI trend f1 and the component D3 is displayed in Fig. 8h. The annual seasonality of the number of fish turns out to lag behind the SOI by about one month, indicating that an increased water temperature induces an increased number of fish about one month later. This effect interferes with the El Niño effect which has periods of abnormal warming of the sea every four to five years. These results confirm similar findings on the interplay between water temperature and fish recruitment by a number of authors such as [8,9,4,10], among others. 6. Final remarks Analyzing multivariate dependence using the nonparametric regression spectrum is particularly useful when the observed series have strong deterministic components. In this article we defined a simple estimator of the regression spectrum based on the periodogram. In contrast to the method in [3] no trend estimation is required. Also, in contrast to the stationary case, no smoothing of the periodogram is needed. In addition, lag estimation in the time domain was considered in order to be able to deal with cases where dependence between two series occurs for a small number of frequencies only. The regression spectrum approach can be particularly powerful when

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

217

Fig. 8. Wavelet estimates of trend components of SOI and the recruitment of new fish are displayed in (a, b). The series in 8b is then decomposed into a long-term component (c) given by the three coarsest levels of wavelet decomposition and a 12-month cycle corresponding to the fourth resolution level (d). The amplitude of cross-spectrum of D4-S6 with SOI is represented in (e) with the solid line given by three times the standard deviation. Important frequencies are again marked by black squares. The estimated phase line in (f) is based on these frequencies. The level D3 essentially corresponds to a 12-month cycle which is represented by a dominating frequency at j = 38 (e). Estimation of the lag by periodogram analysis is inappropriate based on a single frequency. Therefore, the cross-correlations of the D3-component are displayed in (f).

used in combination with multiresolution analysis. Often, the strength, type and interpretation of dependencies differ at different levels. The SOI/fish recruitment data is a typical example of multilevel dependence. In future research, more formal methods should be developed for combining regression spectrum estimation and wavelet decomposition. Acknowledgements This research was supported in part by a grant of the German Research Foundation (DFG Research Unit 512 ‘Price-, Liquidity- and Credit Risks: Measurement and Allocation’). We would like to thank the referees for their very helpful constructive comments and suggestions. Appendix. Proofs Proof 1 (of Theorem 1). We have hrs (j) = ar (j)as (j), and 1 n2

Ar (ωj )As (ωj ) − hrs (j) = O(n−1 ).

(16)

Then 1 n

Irs (ωj ) − hrs (j) =

1 h n2

i

Br (ωj )As (ωj ) + Ar (ωj )Bs (ωj ) + Op (n−1 ) + O(n−1 )

= Op (n−1/2 ).

218

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

For the second part consider E (n−1 Irs (ωj )) =

i 1 h A (ω ) A (ω ) + E ( B (ω ) B (ω )) , r j s j r j s j n2


= hrs (j) + O(n−1 ) + n−1 E (I;rs (ωj )). Results from traditional spectral analysis show that E (I;rs (ωj )) converges to 2π h;rs (ωj ) uniformly for all frequencies ωj so that (2) follows. Note furthermore that 1 

√

n



ar (j)Bs (ωj ) + Br (ωj )as (j) =: α(j) + iβ(j),

(17)

where 1 X

α(j) = √

n

1 X

β(j) = √

n

h

i

h

i

cos(ωj t ) ar (j)s (t ) + as (j)r (t ) ,

(18)

sin(ωj t ) ar (j)s (t ) − as (j)r (t ) .

(19)

Consider now the variance of the adjusted periodogram (the result for covariances follows similarly): var(n−1 Irs (ωj )) = n−1 E [|α(j) + iβ(j)|2 ]

(20)

+ cov(n I;rs (ωj ), α(j) + iβ(j))

(21)

+ cov(α(j) + iβ(j), n I;rs (ωj ))

(22)

+ var(n I;rs (ωj )).

(23)

−1

−1

−1

Using standard results from [2, p. 429] yields

 Σrr Σss + O(n−1 ); cov(I;rs (ωj ), I;rs (ωj0 )) = Σrr Σss + Σrs Σsr + O(n−1 );  −1 O(n );

0 < ωj = ωj0 < π , ωj = ωj0 ∈ {0, π},

ωj 6= ωj0 ,

where the remainders contain the fourth order cumulants between r (i) and s (i). The covariances in (21) and (22) consist of terms of the form 1X

cov(I;rs (ωj ), Br (ωj )) =

n t ,u,v

exp(−iωj (t − u + v))E (r (t )s (u)r (v))

n E (r (1)2 s (1)) X

=

n

exp(iωj t ) = 0,

t =1

where E (r (t )2 s (t )) is independent of t. Simple considerations show that the 2p-dimensional real valued random vector Un (ωj ) = n

X

! (t ) cos(ωj t )

X

(t ) sin(ωj t )

−1/2

(24)

is asymptotically normal with mean 0 and covariance matrix E (Un (ωj )Un (ωj )T ) =

1



2

6 0

0

6



.

Furthermore, for all Fourier frequencies ωj 6= ωj0 , E (Un (ωj )Un (ωj0 )T ) = 0.

(25)

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

219

Therefore,





n−1 E Br (ωj )Bs (ωj )

X = n−1 cov r (t )(cos(ωj t ) − i sin(ωj t )),  X s (t )(cos(ωj t ) − i sin(ωj t )) =

1 2

(Σrs + Σsr ) = Σrs .

Similarly, for all pairs (r , s), 1 ≤ r 6= s ≤ p, var(n−1/2 Br (ωj )) = Σrr

and

cov(Br (ωj ), Bs (ωj )) = cov(Bs (ωj ), Br (ωj )) = 0.

Noting that n X

cos2 (ωj t ) =

t =1

n 2

,

the variance in (20) follows from

τα2 (j) := var(α(j)) n h 1X = cos2 (ωj t ) Σrr |as (j)|2 + Σss |ar (j)|2 + as (j)ar (j)cov(r (t ), s (t )) n t =1

+ ar (j)as (j)cov(s (t ), r (t )) =

1 2

i


Σrr |as (j)|2 + Σss |ar (j)|2 + Σrs Re{as (j)ar (j)}.

Similarly,

τβ2 (j) := var(β(j)) =

1 2


Σrr |as (j)|2 + Σss |ar (j)|2 − Σrs Re{as (j)ar (j)}

and

τα,β (j) := cov(α(j), β(j)) n 1X = cos(ωj t ) sin(ωj t )cov(as (j)r (t ) + ar (j)s (t ), ar (j)s (t ) − as (j)r (t )), n t =1

where the covariance is independent of t. Using the orthogonality relations for trigonometric functions we get that τα,β (j) = 0. For the asymptotic distribution, write

√ n



1 n

 I (ωj ) − hrs (j) = α(j) + iβ(j) + Op (n−1/2 ).

According to (24) the asymptotic distribution of the real and the imaginary part are both univariate normal. Moreover, α(j) and β(j) are uncorrelated and hence asymptotically independent. Hence,

   2 τ (j) α(j) d → N 0, α β(j) 0



0

τβ2 (j)



,

and α(j) + iβ(j) converges in distribution to a complex valued normal random variable with mean 0 and variance 2 2 2 2 2 τα+ iβ (j) = τα (j) + τβ (j) = Σrr |as (j)| + Σss |ar (j)| .

(26)

220

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

This implies

√



n

1 n

 d Irs (ωj ) − hrs (j) → Nc (0, Σrr |as (j)|2 + Σss |ar (j)|2 ).

(27)

The result for finite samples (ωj1 , . . . , ωjk ) follows as usual by applying the Cramer–Wold device.



Proof 2 (of Theorem 2). The following lemma summarizes asymptotic results from [11]. Lemma 2. Assume that the sequence {(i)}, i = 1, . . . , n, is of the form

(i) =

∞ X

A(j)Z(i − j),

(28)

j=−∞

where A(j) = (Alk (j))1≤l,k≤p are p × p-matrices such that for all pairs (l, k), 1

X

|Alk (j)||j| 2 < ∞,

j∈Z

and the sequence Z(i) is independent and identically distributed with mean 0 and non-singular covariance matrix Σ . Then,

X g(ωj ) = n

! (t ) cos(ωj t )

−1/2

X

(t ) sin(ωj t )

g1 (ωj ) g2 (ωj )

 =



with gl (ωj ) = (gl1 (ωj ), . . . , glp (ωj ))T , l ∈ {1, 2}, converges in distribution to a 2p-dimensional random variable with mean 0 and asymptotic covariance matrix

π



C(ωj ) −Q(ωj )

where h (ωj ) =

Q(ωj ) , C(ωj )



    (C(ωj ) − iQ(ωj )) with C(ωj ) = c;rs (ωj ) 1≤r ,s≤p and Q(ωj ) = q;rs (ωj ) 1≤r ,s≤p is the spectral density matrix of (i) at Fourier frequency ωj = 2π j/n. Furthermore g (ωj ) and g (ωj0 ) are asymptotically independent for ωj 6= ωj0 . 1 2

We now turn to the proof of Theorem 2. The first two parts of the proof are similar to those of Theorem 1. For the asymptotic distribution and variance consider Eqs. (20)–(23). Refer again to [2, p. 431] to get cov(I;rs (ωj ), I;rs (ωj0 ))

 (2π )2 h;rr (ωj )h;ss (ωj ) + O(n−1/2 ); = (2π )2 (h;rr (ωj )h;ss (ωj ) + h;rs (ωj )h;sr (ωj )) + O(n−1/2 );  −1 O(n );

0 < ωj = ωj0 < π , ωj = ωj0 ∈ {0, π},

ωj 6= ωj0 .

Using the results of Lemma 2,

 var

 1 √ Br (ωj ) = 2π h;rr (ωj ) (1 ≤ r ≤ p), n

and using the notation h;rr (ωj ) = c;rs (ωj ) − iq;rs (ωj ).

(29)

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

221

Then,

 cov

 1 1 √ Br (ωj ), √ Bs (ωj ) = cov(g1r (ωj ) − ig2r (ωj ), g1s (ωj ) + ig2s (ωj )) n

n

= c;rs (ωj ) + iq;rs (ωj ) − iq;rs (ωj ) − c;rs (ωj ) = 0.

(30)

Due to Lemma 2 and Eqs. (29) and (30), d

α(j) + iβ(j) → Nc (0, (2π |as (j)|2 h;rr (ωj ) + 2π |ar (j)|2 h;ss (ωj ))), where α(j) + iβ(j) is defined in (17). The remaining parts (21) and (22) consist of terms of the form cov(I;rs (ωj ), Br (ωj )) =

P t ,u,v

exp(−iωj (t − u + v))E (r (t )s (u)r (v)),

with p X

E (r (t )s (u)r (v)) =

E (Zr1 (1)Zr2 (1)Zr3 (1))

r1 ,r2 ,r3 =1

X

Ar ,r1 (j)As,r2 (j + u − t )Ar ,r3 (j + v − t ),

j

where Arr1 (j) is the (r , r1 )th entry of the matrix A(j) and E (Zr1 (t )Zr2 (t )Zr3 (t )) is the third order cumulant between the components Zr1 (t ), Zr2 (t ) and Zr3 (t ) of Z(t ). Since (i) is of the form (28),

|cov(I;rs (ωj ), Br (ωj ))| ≤ so that τn (ωj ) = O(n−1 ).

X

|E (r (t )s (u)r (v))| < ∞,

t ,u,v



Proof 3 (of Theorem 3). 1. Results of [6] prove wavelet thresholding estimators to be consistent of order Op (n−1/2 ). The regression cross-covariance is a continuous function of θ , the vector of wavelet coefficient of fr and fs , so that the estimator uˆ max is a continuous function of θˆ . Consistency of θˆ then implies consistency rs of umax . rs 2. Due to (14) and (15),

 
ˆ = 0 and γrs0 umax γrs0 uˆ max rs , θ rs , θ0 = 0. Applying the mean value theorem gives

     
0 max ˆ ˆ , − umax γrs0 uˆ max = γrs00 u˜ rs , θˆ uˆ max rs rs , θ − γrs urs , θ rs and umax where u˜ lies on a line connecting uˆ max rs . Therefore, rs uˆ max − umax =− rs rs

ˆ γrs0 (umax rs ,θ) ˆ . γrs00 (˜urs ,θ)

Furthermore, 0 max ˆ γrs0 (umax rs , θ ) = γrs (urs , θ0 ) + | {z }


T ∂ 0 γ (umax , θ ) θ =θ˜ ∂θ rs rs

· (θˆ − θ0 )

=0

for suitably chosen θ˜ , and uˆ max − umax =− rs rs

1  γrs00 u˜ rs ,θˆ

∂ 0 γ ∂θ rs

umax rs , θ θ=θ˜



T  θˆ − θ0 .

Due to consistency of θˆ and uˆ max rs ,

     00 00 max ˆ + (γrs00 umax ˆ γrs00 (˜urs , θˆ ) = γrs00 (umax urs , θˆ ) − γrs00 umax rs , θ0 ) + γrs (˜ rs , θ rs , θ − γrs (urs , θ0 )) −1/2 ), = γrs00 (umax rs , θ0 ) + Op (n

222

J. Beran, M.A. Heiler / Statistical Methodology 6 (2009) 202–222

and

∂ 0 max ∂ 0 max γrs (urs , θ ) γrs (urs , θ ) + Op (n−1/2 ). = ∂θ ∂θ θ=θ˜ θ=θ0 Then,

T ∂ 0 max ˆ − = − 00 max γ (u , θ0 ) (θˆ − θ0 )(1 + Op (n−1/2 )). γrs (urs , θ0 ) ∂θ rs rs √ √ The central limit theorem for n(θˆ − θ0 ) implies that n(ˆumax − umax rs rs ) is asymptotically normal umax rs



1

umax rs

with mean 0 and covariance

 var(ˆumax rs ) =

∂ 0 γ (umax , θ ) θ=θ ∂θ rs rs 0



T

var(θˆ )



∂ 0 γ (umax , θ ) θ=θ ∂θ rs rs 0

2 γrs00 (umax rs , θ0 )



 . 

References [1] M.B. Priestley, Spectral Analysis and Time Series, 6th ed., Academic Press, London, 1989. [2] P.J. Brockwell, R.A. Davis, Time Series: Theory and Methods, Springer, New York, 1987. [3] J. Beran, M. Heiler, A nonparametric regression spectrum for analyzing multivariate time series, Journal of Multivariate Analysis 99 (4) (2007) 684–714. [4] R.H. Shumway, D.S. Stoffer, Time Series Analysis and its Applications, 2nd ed., Springer, New York, 2000. [5] D. Brillinger, Uses of cumulants in wavelet analysis, Proceedings of SPIE, Advanced Signal Processing 2296 (1994) 2–18. [6] D.R. Brillinger, Some uses of cumulants in wavelet analysis, Nonparametric Statistics 6 (1996) 93–114. [7] D.L. Donoho, I.M. Johnstone, Adapting to unknown smoothness via wavelet shrinkage, Journal of the American Statistical Association 90 (1995). [8] S.A. Murawski, Climate change and marine fish distributions: Forecasting from analogy, Transactions of the American Fisheries Society 122 (1993) 647–658. [9] B.C. Victor, G.M. Wellington, D.R. Robertson, B.I. Ruttenberg, The effect of the EL Niño-Southern Oscillation event on the distribution of reef-associated labrid fishes in the eastern Pacific Ocean, Bulletin of Marine Science 69 (1) (2001) 280–288. [10] O. Rosen, D.S. Stoffer, Automatic estimation of multivariate spectra via smoothing splines, Biometrika 94 (2) (2007) 335–345. [11] E.J. Hannan, Multiple Time Series, Wiley, New York, 1970. [12] J. Beran, S. Ghosh, Estimation of the dominating frequency for stationary and nonstationary fractional autoregressive models, Journal of Time Series Analysis 21 (5) (2000) 517–533.