Chapter 8: Time Series Analysis and Forecasting

CPFFT

TextonSpedometerHPClogo.gif



Computes the cross periodogram of two stationary time series using a fast Fourier transform.

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)

CPM — (N/2 + 1) by 10 matrix containing a summarization of the results 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).

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 series, and the real and imaginary components of the cross periodogram.

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

IPVER — Option for version of the periodogram.   (Input)
Default: IPVER = 0.

IPVER

Action

0

Compute usual periodogram.

1

Compute modified periodogram.

            Refer to the “Description” section for further details.

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).

FORTRAN 90 Interface

Generic:                              CALL CPFFT (X, Y, CPM [,…])

Specific:                             The specific interface names are S_CPFFT and D_CPFFT.

FORTRAN 77 Interface

Single:            CALL CPFFT (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, IPVER, CPM, LDCPM)

Double:                              The double precision name is DCPFFT.

Description

Routine CPFFT computes the cross periodogram of two jointly stationary time series 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 periodogram of the sample sequence {Xt}, t = 1, …, n computed with the padded sequence

is defined by

where

and

represent the

cosine and sine transforms, respectively, and K is the scale factor

The periodogram of the sample sequence {Yt}, t = 1, …, n computed with the padded sequence

is defined 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 periodograms of both

according to the version specified by the argument IPVER. The computational formula for the cross periodogram is given by

where

and

The real part of the (modified) cross periodogram represents the ’raw’ sample cospectrum and the negative of the imaginary part of the (modified) cross periodogram represents the ‘raw’ sample quadrature spectrum (Priestley 1981, page 695). The relationship between the cross periodogram and its complex conjugate is given by

and may be used to recover the cross periodogram at negative frequencies.

Comments

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

CALL C2FFT (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, IPVER, CPM, LDCPM, CX, COEF, WFFTC, CPY)

The additional arguments are as follows:

CX — Complex work vector of length N.

COEF — Complex work vector of length N.

WFFTC — Work vector of length 4N + 15.

CPY — Work vector of length 2N.

2.         The centered and padded time series are defined by

I

AOV(I)

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 cross periodogram IXY(ω) is complex valued in general. The relation
IXY(ω) = conj(IXY(ω)) for w > 0.0 recovers the cross periodogram for negative frequencies since real(IXY(ω)) = real(IXY(ω)) and
imag(IXY(ω)) = imag(IXY(ω)). The periodogram I(ω) is an even function of the frequency ω. The relation I(ω) = I(ω) for ω > 0.0 recovers the periodogram for negative frequencies.

4.         Since cos(ω) is an even function of ω and sin(ω) is an odd function of ω, the cosine and sine transformations, respectively, satisfy A(ω) = A(ω) and B(ω) = B(ω) for
ω > 0.0. Similarly, the complex Fourier coefficients, stored in COEF, satisfy
COEF(ω) = conj(COEF(ω)).

5.         Computation of the 2 * NOBS 1 cross-covariances of X and Y using the inverse Fourier transform of the cross periodogram requires NPAD = NOBS 1.

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 CPFFT to these data produces the following results. Note that CPFFT sets CPM (1, 2) to the missing value code via routine AMACH.p<.STREF.DOC!MACH;;. The printing of CPM (1, 2) depends on the computer.

 

      USE GDATA_INT

      USE CPFFT_INT

      USE WRRRL_INT

 

      IMPLICIT   NONE

      INTEGER    LDCPM, LDRDAT, NDRDAT, NOBS, NPAD

      PARAMETER  (LDRDAT=100, NDRDAT=4, NOBS=100, NPAD=NOBS-1, &

                 LDCPM=(NOBS+NPAD)/2+1)

!

      INTEGER    IPVER, NRCOL, NRROW

      REAL       CPM(LDCPM,10), FLOAT, RDATA(LDRDAT,NDRDAT), &

                 X(NOBS), Y(NOBS)

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

                  TITLE*41

      INTRINSIC  FLOAT

!

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

!

      DATA TITLE/'Results of the Cross Periodogram Analysis'/

      DATA FMT/'(F10.3)'/

      DATA CLABEL1/'k+1', 'w(k)', 'p(k)', 'IX(w(k))', 'AX(w(k))', &

          'BX(w(k))'/

      DATA CLABEL2/'k+1', 'IY(w(k))', 'AY(w(k))', 'BY(w(k))', &

          'Real IXY', 'Imag. IXY'/

      DATA RLABEL/'NUMBER'/

!

!                                 Robinson Data

      CALL GDATA (8, RDATA, NRROW, NRCOL)

!                                 Center on arithmetic means

!                                 Frequency in radians per unit time

!                                 Modified periodogram version

      IPVER = 1

!                                 Compute the cross periodogram

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

!

!                                 Print results (First 10 rows)

      CALL WRRRL (TITLE, CPM, RLABEL, CLABEL1, 10, 5, FMT=FMT)

      CALL WRRRL ('%/', CPM(1:,6), RLABEL, CLABEL2, 10, 5, FMT=FMT)

!

      END

Output

 

          Results of the Cross Periodogram Analysis
k+1        w(k)        p(k)    IX(w(k))    AX(w(k))    BX(w(k))
  1       0.000         NaN       0.000       0.000       0.000
  2       0.032     199.000     184.159       3.742     -13.044
  3       0.063      99.500    1364.408      35.457     -10.354
  4       0.095      66.333    2433.933      29.411      39.610
  5       0.126      49.750    1351.002     -21.749      29.631
  6       0.158      39.800     140.421     -11.716      -1.773
  7       0.189      33.167      44.117      -4.671       4.722
  8       0.221      28.429     121.186     -11.003      -0.343
  9       0.253      24.875     176.275      -4.782     -12.386
 10       0.284      22.111     144.867      10.038      -6.642

k+1    IY(w(k))    AY(w(k))    BY(w(k))    Real IXY   Imag. IXY
  1       0.000       0.000       0.000       0.000       0.000
  2    1689.212     -37.480     -16.866      79.776    -552.014
  3    4113.003      41.232     -49.122    1970.577   -1314.779
  4    3255.785      44.214      36.068    2729.031    -690.474
  5    1757.663      -8.162      41.122    1396.006    -652.513
  6    1002.050     -30.107       9.778     335.410    -167.954
  7      62.360      -6.825       3.972      50.636      13.678
  8    1481.396     -38.096       5.487     417.288     -73.451
  9    1274.161     -17.176     -31.291     469.704     -63.095
 10     488.479     -12.442     -18.267      -3.570    -265.992



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