Chapter 8: Time Series Analysis and Forecasting

CSSWD

TextonSpedometerHPClogo.gif



Estimates the nonnormalized cross-spectral density of two stationary time series using a spectral window given the time series data.

Required Arguments

X — Vector of length NOBS containing the first stationary time series.   (Input)

Y — Vector of length NOBS containing the second stationary time series.   (Input)

F — Vector of length NF containing the frequencies at which to evaluate the cross-spectral density estimate.   (Input)
The units of F correspond to the scale specified by IFSCAL. The elements of F must be in the range (π/TINT, π /TINT), inclusive, for IFSCAL = 0 and
(1/(2 * TINT), 1/(2 * TINT)), inclusive, for IFSCAL = 1.

ISWVER — Option for version of the spectral window.   (Input)

ISWVER

Action

1

Modified Bartlett

2

Daniell

3

Tukey-Hamming

4

Tukey-Hanning

5

Parzen

6

Bartlett-Priestley

            Refer to the “Algorithm” section for further details.

M — Vector of length NM containing the values of the spectral window parameter M.   (Input)
For the Parzen spectral window (ISWVER = 5), all values of the spectral window parameters M must be even.

CPM — (N/2 + 1) by 10 matrix containing a summarization of the cross periodogram analysis.   (Output)
For k = 0, 1, …, N/2, the (k + 1)-st element of the j-th column of CPM is defined as

Col.

Description

1

Frequency, ωk where ωk = 2πk/N for IFSCAL = 0 or ωk = k/N for IFSCAL = 1.

2

Period, pk where pk = 2π/ωk for IFSCAL = 0 and pk = 1/ωk for IFSCAL = 1. If ωk = 0, pk is set to missing.

3

X periodogram ordinate, IXk)

4

X cosine transformation coefficient, AXk)

5

X sine transformation coefficient, BXk)

6

Y periodogram ordinate, IYk)

7

Y cosine transformation coefficient, AYk)

8

Y sine transformation coefficient, BYk)

9

Real part of the XY cross periodogram ordinate IXYk).

10

Imaginary part of the XY cross periodogram ordinate IXYk).

            Note N = NOBS + NPAD.

CSM — NF by (NM * 7 + 2) matrix containing a summarization of the cross-spectral analysis.   (Output)
The k-th element of the j-th column of CSM is defined as

Col.

Description

1

Frequency, F(k).

2

Period, pk where pk = 2π/F(k) for IFSCAL = 0 and pk = 1/F(k) for IFSCAL = 1. If F(k) = 0, pk is set to missing.

3

spectral density estimate at F(k) using the spectral window parameter M(1).

4

Y spectral density estimate at F(k) using the spectral window parameter M(1).

5

Cospectrum estimate at F(k) using the spectral window parameter M(1).

6

Quadrature spectrum estimate at F(k) using the spectral window parameter M(1).

7

Cross-amplitude spectrum estimate at F(k) using the spectral window parameter M(1).

8

Phase spectrum estimate at F(k) using the spectral window parameter M(1).

9

Coherence estimate at F(k) using the spectral window parameter M(1).

 

NM * 7 + 2

Coherence estimate at F(k) using the spectral window parameter M(NM).

            where k = 1, …, NF.

Optional Arguments

NOBS — Number of observations in each stationary time series X and Y.   (Input)
NOBS must be greater than or equal to two.
Default: NOBS = size (X,1).

IPRINT — Printing option.   (Input)
Default: IPRINT = 0.

IPRINT

Action

0

No printing is performed.

1

Prints the cross periodogram and cross-spectral density estimate based on a specified version of a spectral window for a given set of spectral window parameters.

XCNTR — Constant used to center the time series X.   (Input)
Default: XCNTR = the arithmetic mean.

YCNTR — Constant used to center the time series Y.   (Input)
Default: YCNTR = the arithmetic mean.

NPAD — Number of zeroes used to pad each centered time series.   (Input)
NPAD must be greater than or equal to zero. The length of each centered and padded time series is N = NOBS + NPAD.
Default: NPAD = NOBS – 1.

IFSCAL — Option for frequency scale.   (Input)
Default: IFSCAL = 0.

IFSCAL

Action

0

Frequency in radians per unit time.

1

Frequency in cycles per unit time.

NF — Number of frequencies at which to evaluate the cross-spectral density estimate.   (Input)
Default: NF = size (F,1).

TINT — Time interval at which the series are sampled.   (Input)
For a discrete parameter process, usually TINT = 1. For a continuous parameter process, TINT > 0. TINT is used to adjust the cross-spectral density estimate.
Default: TINT = 1.0.

NM — Number of spectral window parameters M used to compute the cross-spectral density estimate for a given spectral window version.   (Input)
NM must be greater than or equal to one.
Default: NM = size (M,1).

LDCPM — Leading dimension of CPM exactly as specified in the dimension statement of the calling program.   (Input)
LDCPM must be greater than or equal to N/2, + 1.
Default: LDCPM = size (CPM,1).

LDCSM — Leading dimension of CSM exactly as specified in the dimension statement of the calling program.   (Input)
LDCSM must be greater than or equal to NF.
Default: LDCSM = size (CSM,1).

FORTRAN 90 Interface

Generic:                              CALL CSSWD (X, Y, F, ISWVER, M, CPM, CSM[,…])

Specific:                             The specific interface names are S_CSSWD and D_CSSWD.

FORTRAN 77 Interface

Single:            CALL CSSWD (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, ISWVER, NM, M, CPM, LDCPM, CSM, LDCSM)

Double:                              The double precision name is DCSSWD.

Description

Routine CSSWD estimates the nonnormalized cross-spectral density function of two jointly stationary time series using a spectral window given a sample of n = NOBS observations {Xt} and {Yt} for t = 1, 2, …, n.

Let

represent the centered and padded data where N = NOBS + NPAD,

and

is determined by

Similarly, let

represent the centered and padded data where

and

is determined by

The modified periodogram of

is estimated by

where

and

represent the

cosine and sine transforms, respectively, and K is the scale factor equal to 1/(2πn). The modified periodogram of

is estimated by

where

and

represent the

cosine and sine transforms, respectively. Since the periodogram is an even function of the frequency, it is sufficient to estimate the periodogram at the discrete set of nonnegative frequencies

The routine PFFT is used to compute the modified periodograms of both

The computational formula for the cross periodogram is given by

where

and

The routine CPFFT is used to compute the modified cross periodogram between

The nonnormalized spectral density of Xt is estimated by

and the nonnormalized spectral density of Yt is estimated by

where the spectral window Wn(θ) is specified by argument ISWVER. The following spectral windows Wn(θ) are available.

Modified Bartlett

where FM(θ) corresponds to the Fejιr kernel of order M.

Daniell

Tukey

where DM(θ) represents the Dirichlet kernel. The Tukey-Hamming window is obtained when
a = 0.23, and the Tukey-Hanning window is obtained when a = 0.25.

Parzen

where M is even. If M is odd, then M + 1 is used instead of M in the above formula.

Bartlett-Priestley

The argument NM specifies the number of window parameters M and, hence, corresponds to the number of spectral density estimates to be computed for a given spectral window. Note that the same spectral window Wn(θ) and set of parameters M are used to obtain both

The above spectral density formulas assume the data {Xt} and {Yt} correspond to a realization of a bivariate discrete-parameter stationary process observed consecutively in time. In this case, the observations are equally spaced in time with interval Δt = TINT equal to one. However, if the data correspond to a realization of a bivariate continuous-parameter stationary process recorded at equal time intervals, then the spectral density estimates must be adjusted for the effect of aliasing. In general, the estimate of hX(ω) is given by

and the estimate of hY(ω) is given by

The nonnormalized spectral density is estimated over the set of frequencies

ω = fi,     i = 1, …, nf

where nf = NF. These frequencies are in the scale specified by the argument IFSCAL but are transformed to the scale of radians per unit time for computational purposes. The frequency ω of the desired spectral estimate is assumed to be input in a form already adjusted for the time interval Δt.

The cross-spectral density function is complex-valued in general and may be written in the following form:

The cospectrum is estimated by

and the quadrature spectrum is estimated by

Note that the same spectral window Wn(θ) and window parameter M used to derive

are also used to compute

The nonnormalized cross-spectral density estimate is computed over the same set of frequencies as the nonnormalized spectral density estimates with a similar adjustment for Δt.

An equivalent representation of hXY(ω) is the polar form defined by

The cross-amplitude spectrum is estimated by

and the phase spectrum is estimated by

Finally, the coherency spectrum is estimated by

The coherence or squared coherency is output.

Comments

1.         Workspace may be explicitly provided, if desired, by use of C2SWD/DC2SWD. The reference is:

CALL C2SWD (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, ISWVER, NM, M, CPM, LDCPM, CSM, LDCSM, CX, COEF, WFFTC, CPY)

The additional arguments are as follows:

CX — Complex work vector of length N.   (Output)

COEF — Complex work vector of length N.   (Output)

WFFTC — Vector of length 4N + 15.

CPY — Vector of length 2N.

2.         The centered and padded time series are defined by

CX(j) = X(j)XCNTR for j = 1, …, NOBS

CX(j) = 0       for j = NOBS + 1, …, N

and

CY(j) = Y(j)YCNTR  for j = 1, …, NOBS

CY(j) = 0       for j = NOBS + 1, …, N

where N = NOBS + NPAD.

3.         The normalized cross-spectral density estimate is obtained by dividing the nonnormalized cross-spectral density estimate in matrix CSM by the product of the estimated standard deviation of X and the estimated standard deviation of Y.

Example

Consider the Robinson Multichannel Time Series Data (Robinson 1967, page 204) where X is the Wφlfer sunspot number and Y is the northern light activity for the time period from 1770 through 1869. Application of routine CSSWD to these data produces the following results:

 

      USE UMACH_INT

      USE GDATA_INT

      USE CSSWD_INT

      USE WRRRL_INT

 

      IMPLICIT   NONE

      INTEGER    LDCPM, LDCSM, LDRDAT, N, NDRDAT, NF, NM, &

                NOBS, NPAD

      PARAMETER  (LDRDAT=100, NDRDAT=4, NF=10, NM=2, &

                  NOBS=100, LDCSM=NF, NPAD=NOBS-1, N=NOBS+NPAD,&

                  LDCPM=N/2+1)

!

      INTEGER    I, ISWVER, J, JPT, M(NM), NOUT, NRCOL, NRROW

      REAL       ASIN, CPM(LDCPM,10), CSM(LDCSM,NM*7+2), F(NF), FLOAT, &

                 PI, RDATA(LDRDAT,NDRDAT), TINT, X(NOBS), Y(NOBS)

      CHARACTER  CLABEL1(3)*9, CLABEL2(6)*16, FMT*7, RLABEL(1)*6, &

                 TITLE*80

      INTRINSIC  ASIN, FLOAT

!

      EQUIVALENCE (X(1), RDATA(1,2)), (Y(1), RDATA(1,3))

!

      DATA FMT/'(F10.4)'/

      DATA CLABEL1/' k', 'Frequency', 'Period'/

      DATA CLABEL2/'%/ k', '%/Cospectrum', '%/Quadrature', &

          'Cross%/Amplitude', '%/Phase', '%/Coherence'/

      DATA RLABEL/'NUMBER'/

!                                 Initialization

      CALL UMACH (2, NOUT)

      PI = 2.0*ASIN(1.0)

      DO 10  I=1, NF

         F(I) = PI*FLOAT(I)/FLOAT(NF)

   10 CONTINUE

!                                 Robinson Data

      CALL GDATA (8, RDATA, NRROW, NRCOL)

!                                 Center on arithmetic means

!                                 Frequency in radians per unit time

!                                 Spectral window parameters

      M(1) = 10

      M(2) = 30

!                                 Time interval for discrete data

!                                 Compute cross-spectral density

!                                 using the Parzen window

      ISWVER = 5

      CALL CSSWD (X, Y, F, ISWVER, M, CPM, CSM)

!                                 Print results

      TITLE = 'Cross-Spectral Analysis Using Parzen Window'

      CALL WRRRL (TITLE, CSM, RLABEL, CLABEL1, NF, 2, FMT=FMT)

      DO 20  J=1, NM

         JPT   = 7*(J-1) + 5

         TITLE = '%/Results of the Cross-Spectral Analysis With '// &

                'Spectral Window Parameter M = '

         WRITE (TITLE(77:78),'(I2)') M(J)

         CALL WRRRL (TITLE, CSM(1:,JPT:), RLABEL, CLABEL2, NF, 5, FMT=FMT)

   20 CONTINUE

!

      END

Output

 

Cross-Spectral Analysis Using Parzen Window
         k   Frequency      Period
         1      0.3142     20.0000
         2      0.6283     10.0000
         3      0.9425      6.6667
         4      1.2566      5.0000
         5      1.5708      4.0000
         6      1.8850      3.3333
         7      2.1991      2.8571
         8      2.5133      2.5000
         9      2.8274      2.2222
        10      3.1416      2.0000

Results of the Cross-Spectral Analysis With Spectral Window Parameter M = 10
                                         Cross
         k  Cospectrum  Quadrature   Amplitude       Phase   Coherence
         1    463.5888    -65.9763    468.2600      0.1414      0.2570
         2    286.5450    -75.0209    296.2029      0.2561      0.1710
         3    150.1073    -57.8263    160.8604      0.3677      0.1438
         4     52.9840    -32.3642     62.0866      0.5483      0.0998
         5     21.5435    -15.0888     26.3020      0.6110      0.0794
         6     21.4228     -9.8188     23.5658      0.4298      0.1716
         7     15.7005     -5.3704     16.5936      0.3296      0.2112
         8      8.0118     -1.8887      8.2314      0.2315      0.1272
         9      2.7682      0.2007      2.7754     -0.0724      0.0446
        10      0.5777      0.1008      0.5864     -0.1727      0.0091

Results of the Cross-Spectral Analysis With Spectral Window Parameter M = 30
                                         Cross
         k  Cospectrum  Quadrature   Amplitude       Phase   Coherence
         1    169.7542   -193.4384    257.3615      0.8505      0.1620
         2    452.6187     32.3813    453.7755     -0.0714      0.2213
         3     94.5221    -90.8159    131.0800      0.7654      0.2629
         4     -0.2096     -6.1127      6.1163      1.6051      0.0019
         5     27.4711    -22.1946     35.3166      0.6796      0.2492
         6     29.1329     -4.0128     29.4080      0.1369      0.3170
         7     11.2058     -9.3403     14.5881      0.6948      0.2594
         8      8.0017      0.8813      8.0501     -0.1097      0.1928
         9     -0.4199      2.2893      2.3275     -1.7522      0.0468
        10      0.5570     -1.0767      1.2123      1.0934      0.0678



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260