Chapter 8: Time Series Analysis and Forecasting

CSSWP

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

Required Arguments

N — Number of observations in each of the appropriately centered and padded time series X and Y.   (Input)
N must be greater than or equal to two.

SX — Vector of length NF containing the estimate of the spectral density of the first time series X.   (Input)

SY — Vector of length NF containing the estimate of the spectral density of the second time series Y.   (Input)

CPREAL — Vector of length N/2 + 1 containing the real part of the cross periodogram between X and Y.   (Input)
The real part of the cross periodogram evaluated at (angular) frequency wk = 2πk/N is given by CPREAL(k + 1), k = 0, 1, …, N/2.

CPIMAG — Vector of length N/2 + 1 containing the imaginary part of the cross periodogram between X and Y.   (Input)
The imaginary part of the cross periodogram evaluated at (angular) frequency
wk = 2πk/N is given by CPIMAG(k + 1), k = 0, 1, …, N/2.

F — Vector of length NF containing the (angular) frequencies at which the spectral and cross-spectral densities are estimated.   (Input)

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

SWVER

Action

1

Modified Bartlett

2

Daniell

3

Tukey-Hamming

4

Tukey-Hanning

5

Parzen

6

Bartlett-Priestley

            Refer to the “Description” section for further details.

M — Spectral window parameter.   (Input)
M must be greater than or equal to one and less than N. For the Parzen spectral window (ISWVER = 5), the spectral window parameter M must be even.

COSPEC — Vector of length NF containing the estimate of the cospectrum.   (Output)

QUADRA — Vector of length NF containing the estimate of the quadrature spectrum.   (Output)

CRAMPL — Vector of length NF containing the estimate of the cross-amplitude spectrum.   (Output)

PHASE — Vector of length NF containing the estimate of the phase spectrum.   (Output)

COHERE — Vector of length NF containing the estimate of the coherence or squared coherency.   (Output)

Optional Arguments

NF — Number of (angular) frequencies.   (Input)
NF must be greater than or equal to one. Default: NF = size (F,1).

FORTRAN 90 Interface

Generic:          CALL CSSWP (N, SX, SY, CPREAL, CPIMAG, F, ISWVER, M, COSPEC, QUADRA, CRAMPL, PHASE, COHERE [,…])

Specific:                             The specific interface names are S_CSSWP and D_CSSWP.

FORTRAN 77 Interface

Single:            CALL CSSWP (N, SX, SY, CPREAL, CPIMAG, NF, F, ISWVER, M, COSPEC, QUADRA, CRAMPL, PHASE, COHERE)

Double:                              The double precision name is DCSSWP.

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 CSSWP to these data produces the following results.

 

      USE IMSL_LIBRARIES

 

      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, IPVER, ISWVER, J, JPT, JST, M(NM), NRCOL, NRROW

      REAL       COHERE(NF), COSPEC(NF), CPIMAG(LDCPM), &

                 CPM(LDCPM,10), CPREAL(LDCPM), CRAMPL(NF), &

                 CSM(LDCSM,7*NM+2), F(NF), FLOAT, P(NF), PHASE(NF), &

                 PI, PX(LDCPM), PY(LDCPM), QUADRA(NF), &

                 RDATA(LDRDAT,NDRDAT), SX(NF), SY(NF), X(NOBS), Y(NOBS)

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

                 TITLE*80

      INTRINSIC  FLOAT

!

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

      EQUIVALENCE (PX(1), CPM(1,3)), (PY(1), CPM(1,6))

      EQUIVALENCE (CPREAL(1), CPM(1,9)), (CPIMAG(1), CPM(1,10))

      EQUIVALENCE (CSM(1,1), F(1)), (CSM(1,2), P(1))

!

      DATA FMT/'(F12.4)'/

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

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

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

      DATA RLABEL/'NUMBER'/

!                                 Initialization

      PI = 2.0*ASIN(1.0)

      DO 10  I=1, NF

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

         P(I) = 2.0*FLOAT(NF)/FLOAT(I)

   10 CONTINUE

!                                 Robinson Data

      CALL GDATA (8, RDATA, NRROW, NRCOL)

!                                 Center on arithmetic means

!                                 Frequency in radians per unit time

!                                 Modified periodogram version

      IPVER = 1

!                                 Compute cross periodogram

      CALL CPFFT (X, Y, CPM, IPVER=IPVER)

!                                 Spectral window parameters

      M(1) = 10

      M(2) = 30

!                                 Compute cross-spectral density

!                                 using the Parzen window

!

!                                 Print frequency and period

      TITLE = 'Cross-Spectral Analysis Using Parzen Window'

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

      ISWVER = 5

      DO 20  J=1, NM

!                                 Estimate the spectral densities

         CALL SSWP (N, PX, F, M(J), SX, ISWVER=ISWVER)

         CALL SSWP (N, PY, F, M(J), SY, ISWVER=ISWVER)

!                                 Estimate the cross-spectral density

         CALL CSSWP (N, SX, SY, CPREAL, CPIMAG, F, ISWVER, M(J), &

                    COSPEC, QUADRA, CRAMPL, PHASE, COHERE)

!                                 Copy results to output matrices

         JPT = 7*(J-1) + 2

         JST = 7*(J-1) + 5

         CALL SCOPY (NF, SX, 1, CSM(1:,JPT+1), 1)

         CALL SCOPY (NF, SY, 1, CSM(1:,JPT+2), 1)

         CALL SCOPY (NF, COSPEC, 1, CSM(1:,JPT+3), 1)

         CALL SCOPY (NF, QUADRA, 1, CSM(1:,JPT+4), 1)

         CALL SCOPY (NF, CRAMPL, 1, CSM(1:,JPT+5), 1)

         CALL SCOPY (NF, PHASE, 1, CSM(1:,JPT+6), 1)

         CALL SCOPY (NF, COHERE, 1, CSM(1:,JPT+7), 1)

!                                 Print results

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

                'Spectral Window Parameter M = '

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

         CALL WRRRL (TITLE, CSM(1:,JST), 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

Description

Routine CSSWP estimates the nonnormalized cross-spectral density function of two jointly stationary time series using a spectral window given the modified cross-periodogram and spectral densities of the appropriately centered and padded data

for t = 1, …, N.

The routine CPFFT may be used to compute the modified periodograms

and cross periodogram

over the discrete set of nonnegative frequencies

(Here, a means the greatest integer less than or equal to a.) Either routine SSWP or routine SWEP may be applied to the periodograms to obtain nonnormalized spectral density estimates

over the set of frequencies

ω = fi,      i = 1, …, nf

where nf = NF. These frequencies are in the scale of radians per unit time. The time sampling interval Δt is assumed to be equal to one. Note that the spectral window or weight sequence used to compute

may differ from that used to compute

The cross-spectral density function is complex-valued in general and may be written as

The cospectrum is estimated by

and the quadrature spectrum 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

Only one window parameter M may be specified so that only one estimate of hXY (ω) is computed. The nonnormalized cross-spectral density estimate is computed over the same set of frequencies as the nonnormalized spectral density estimates discussed above. However, the particular spectral window used to compute

need not correspond to either the spectral window or the weight sequence used to compute either

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.         The periodograms of X and Y and cross periodogram between X and Y may be computed using the routine CPFFT. The spectral densities of X and Y may then be estimated using any of the routines SSWD, SWED, SSWP, or SWEP. Thus, different window types and/or weight sequences may be used to estimate the spectral and cross-spectral densities given either the series or their periodograms. Note that use of the modified periodograms and modified cross periodogram ensures that the scale of the spectral and cross-spectral densities and their estimates is equivalent.

2          The time sampling interval, TINT, is assumed to be equal to one. This assumption is appropriate for discrete parameter processes. The adjustment for continuous parameter processes (TINT > 0.0) involves multiplication of the frequency vector F by 1/TINT and multiplication of the spectral and cross-spectral density estimates by TINT.

3.         To convert the frequency scale from radians per unit time to cycles per unit time, multiply F by 1/(2π).



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