Chapter 8: Time Series Analysis and Forecasting

CSWED

TextonSpedometerHPClogo.gif



Estimates the nonnormalized cross-spectral density of two stationary time series using a weighted cross periodogram 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.

WT — Vector of length NWT containing the weights used to smooth the periodogram.   (Input)
The actual weights are the values in WT normalized to sum to 1 with the current periodogram ordinate taking the middle weight for NWT odd or the weight to the right of the middle for NWT 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).

CSM — NF by 9 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

X spectral density estimate at F(k) using the specified relative weights contained in WT.

4

Y spectral density estimate at F(k) using the specified relative weights contained in WT.

5

Co-spectrum estimate at F(k) using the specified relative weights contained in WT.

6

Quadrature spectrum estimate at F(k) using the specified relative weights contained in WT.

7

Cross-amplitude spectrum estimate at F(k).

8

Phase spectrum estimate at F(k).

9

Coherence estimate at F(k).

            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 periodogram, cosine and sine transformations of each centered  and padded time series, the real and imaginary components of the cross periodogram, and the cross-spectral density estimate based on a specified weight sequence.

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.0. For a continuous parameter process, TINT > 0.0. TINT is used to adjust the cross-spectral density estimate.
Default: TINT = 1.0.

NWT — Number of weights.   (Input)
NWT must be greater than or equal to one.
Default: NWT= size (WT,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 CSWED (X, Y, F, WT, CPM, CSM [,…])

Specific:                             The specific interface names are S_CSWED and D_CSWED.

FORTRAN 77 Interface

Single:            CALL CSWED (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, NWT, WT, CPM, LDCPM, CSM, LDCSM)

Double:                              The double precision name is DCSWED.

Description

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

for t = 1, …, N represent the centered and padded data where N = NOBS + NPAD,

and

is determined by

Similarly, let

for t = 1, …, N represent the centered and padded data where

and

is determined by

The modified periodogram of

for t = 1, …, N 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 {Yt} for t = 1, …, N 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

(Here, a means the greatest integer less than or equal to a). 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

and k(ω) is the integer such that ωk,0 is closest to ω. The sequence of m = NWT weights {wj} for
j = m/2, …, (m m/2 1) satisfies Σjwj = 1. These weights are fixed in the sense that they do not depend on the frequency ω at which to estimate the spectral density. Usually, m is odd with the weights symmetric about the middle weight w0. If m is even, the weight to the right of the middle is considered w0. The argument WT may contain relative weights since they are normalized to sum to one in the actual computations. 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 equivalent 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

ω = ƒi,     i = 1, …, nƒ

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 as

hXY(ω) = cXY (ω) iqXY(ω)

The cospectrum is estimated by

and the quadrature spectrum is estimated by

Note that the same sequence of weights {wj} used to estimate

is used to estimate

The nonnormalized cross-spectral density estimate is computed over the same set of frequencies as the nonnormalized spectral density estimates discussed above 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 C2WED/DC2WED. The reference is:

CALL C2WED (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, NF, F, TINT, NWT, WT, CPM, LDCPM, CSM, LDCSM, CWK, COEFWK, WFFTC, CPY)

The additional arguments are as follows:

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

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

WFFTC — Vector of length 4N + 15.

CPY — Vector of length 2N.

2.         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 years 1770 through 1869. Application of routine CSWED to these data produces the following results.

 

      USE IMSL_LIBRARIES

 

      IMPLICIT   NONE

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

                 NPAD, NWT

      PARAMETER  (LDRDAT=100, NDRDAT=4, NF=10, NOBS=100, &

                 NWT=7, LDCSM=NF, NPAD=NOBS-1, N=NPAD+NOBS, &

                 LDCPM=N/2+1)

!

      INTEGER    I, NRCOL, NRROW

      REAL       CPM(LDCPM,10), CSM(LDCSM,9), F(NF), FLOAT, PI, &

                 RDATA(LDRDAT,NDRDAT), WT(NWT), X(NOBS), Y(NOBS)

      CHARACTER  CLABEL1(5)*24, CLABEL2(6)*16, FMT*7, RLABEL(1)*6, &

                 TITLE1*32, TITLE2*40

      INTRINSIC  FLOAT

!

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

!

      DATA WT/1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0/

      DATA FMT/'(F12.4)'/

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

          'Spectral%/Estimate%/of X', 'Spectral%/Estimate%/of Y'/

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

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

      DATA RLABEL/'NUMBER'/

      DATA TITLE1/'Results of the Spectral Analyses'/

      DATA TITLE2/'%/Results of the Cross-Spectral Analysis'/

!                                 Initialization

      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

!                                 Time interval for discrete data

!                                 Compute the cross periodogram

      CALL CSWED (X, Y, F, WT, CPM, CSM)

!                                 Print results

      CALL WRRRL (TITLE1, CSM, RLABEL, CLABEL1, NF, 4, FMT=FMT)

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

!

      END

Output

 

             Results of the Spectral Analyses
                                    Spectral      Spectral
                                    Estimate      Estimate
 k     Frequency        Period          of X          of Y
 1        0.3142       20.0000      116.9550     1315.8370
 2        0.6283       10.0000     1206.6086     1005.1219
 3        0.9425        6.6667       84.8369      317.2589
 4        1.2566        5.0000       55.2120      270.2111
 5        1.5708        4.0000       46.5748      115.6768
 6        1.8850        3.3333       12.4050      250.0125
 7        2.1991        2.8571        7.0934       82.6773
 8        2.5133        2.5000        3.4091       62.3267
 9        2.8274        2.2222        5.6828       12.8970
10        3.1416        2.0000        4.0346       17.6441

                 Results of the Cross-Spectral Analysis
                                       Cross
 k    Cospectrum    Quadrature     Amplitude         Phase     Coherence
 1       94.0531     -254.0125      270.8659        1.2162        0.4767
 2      702.5118       21.9823      702.8557       -0.0313        0.4073
 3       70.2379      -31.4431       76.9547        0.4209        0.2200
 4       -1.8715      -36.1639       36.2123        1.6225        0.0879
 5       36.6366      -18.5925       41.0843        0.4696        0.3133
 6       32.7071       -6.6569       33.3776        0.2008        0.3592
 7        9.4887       -9.1692       13.1950        0.7683        0.2969
 8        9.2534       -0.3000        9.2583        0.0324        0.4034
 9       -0.5568        2.9455        2.9977       -1.7576        0.1226
10        1.7640       -1.8321        2.5433        0.8043        0.0909



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