Chapter 8: Time Series Analysis and Forecasting

SSWP

Estimates the nonnormalized spectral density of a stationary time series using a spectral window given the periodogram.

Required Arguments

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

PX — Vector of length N/2 + 1 containing the (modified) periodogram of X.   (Input)
The periodogram ordinate evaluated at (angular) frequency wk = 2πk/N is given by
PX(k + 1), k = 0, 1, …, N/2.

F — Vector of length NF containing the (angular) frequencies at which the spectral density is estimated.   (Input)

M — Spectral window parameter.   (Input)
M must be greater than or equal to one and less than N.

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

Optional Arguments

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

ISWVER — Option for version of the spectral window.   (Input)
Default: ISWVER = 1.

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.

FORTRAN 90 Interface

Generic:                              CALL SSWP (N, PX, F, M, SX [,…])

Specific:                             The specific interface names are S_SSWP and D_SSWP.

FORTRAN 77 Interface

Single:                                CALL SSWP (N, PX, NF, F, ISWVER, M, SX)

Double:                              The double precision name is DSSWP.

Description

Routine SSWP estimates the nonnormalized spectral density function of a stationary time series using a spectral window given the modified periodogram of the appropriately centered and padded data

The routine PFFT may be used to obtain the modified periodogram

over the discrete set of nonnegative frequencies

The symmetry of the periodogram is used to recover the ordinates at negative frequencies.

The estimate of the nonnormalized spectral density hX(ω) is computed according to

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

Modified Bartlett

where FM(θ) corresponds to the Fejer 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 hX(ω) is computed. The nonnormalized spectral density is estimated over the set of frequencies

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

where nƒ = NF. These frequencies are in the scale of radians per unit time. The time sampling interval Δt is assumed to be equal to one.

Comments

1.         The periodogram of X may be computed using the routine PFFT. Estimation of the spectral density of X using the modified periodogram preserves the scale of the spectral density up to adjustment for the time sampling interval.

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 density estimate by TINT.

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

Example

Consider the W๖lfer Sunspot Data (Anderson 1971, page 660) consisting of the number of sunspots observed each year from 1749 through 1924. The data set for this example consists of the number of sunspots observed from 1770 through 1869. Application of routine SSWP to these data produces the following results:

 

      USE IMSL_LIBRARIES

 

      IMPLICIT   NONE

      INTEGER    LDPM, LDSM, NF, NM, NOBS

      REAL       PI, NPAD

      PARAMETER  (NF=20, NM=3, NOBS=100, PI=3.141592654, &

                 LDPM=NOBS, LDSM=NF)

!

      INTEGER    I, IPVER, ISWVER, J, M(NM), N, NCOL, NROW

      REAL       F(NF), PM(LDPM,5), PX(LDPM), RDATA(176,2), FLOAT, &

                 SM(NF,5), SX(NF), X(NOBS)

      CHARACTER  CLABEL(6)*9, FMT*20, RLABEL(1)*6, TITLE*60

      INTRINSIC  FLOAT

!

      EQUIVALENCE (PX(1), PM(1,3)), (F(1), SM(1,1))

      EQUIVALENCE (X(1), RDATA(22,2))

!

      DATA RLABEL/'NONE'/, CLABEL/' ', 'Frequency', 'Period', &

          'M = 10', 'M = 20', 'M = 30'/

!                                 Wolfer Sunspot Data for

!                                 years 1770 through 1869

      CALL GDATA (2, RDATA, NROW, NCOL)

!                                 Center on arithmetic mean

!                                 Pad standard amount

      NPAD = NOBS-1

!                                 Frequency in radians per unit time

!                                 Modified periodogram version

      IPVER = 1

!                                 Compute periodogram

      CALL PFFT (X, PM, IPVER=IPVER)

!                                 Number of observations used to

!                                 compute the periodogram

      N = NOBS + NPAD

!                                 Determine frequency and period

!                                 at which to evaluate the spectral

!                                 density

      DO 10  I=1, NF

         SM(I,1) = PI*FLOAT(I)/FLOAT(NF)

         SM(I,2) = 2.0*FLOAT(NF)/FLOAT(I)

   10 CONTINUE

!                                 Spectral window parameters

      M(1) = 10

      M(2) = 20

      M(3) = 30

!                                 Compute spectral density using

!                                 the Parzen window

      ISWVER = 5

      DO 20  J=1, NM

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

!                                 Copy into SM

         CALL SCOPY (NF, SX, 1, SM(1:,2+J), 1)

   20 CONTINUE

!                                 Print results

      TITLE = 'Spectral Density Using the Parzen Window'

      FMT   = '(F9.4, F6.2, 3F10.2)'

      CALL WRRRL (TITLE, SM, RLABEL, CLABEL, FMT=FMT)

!                                 Compute spectral density using

!                                 the Bartlett-Priestley window

      ISWVER = 6

      DO 30  J=1, NM

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

!                                 Copy into SM

         CALL SCOPY (NF, SX, 1, SM(1:,2+J), 1)

   30 CONTINUE

!                                 Print results

      TITLE = '%/Spectral Density Using the Bartlett-Priestley '// &

             & 'Window'

      CALL WRRRL (TITLE, SM, RLABEL, CLABEL, FMT=FMT)

!

      END

Output

 

      Spectral Density Using the Parzen Window
Frequency  Period      M = 10      M = 20      M = 30
   0.1571   40.00      659.64      617.42      619.73
   0.3142   20.00      666.95      554.70      339.61
   0.4712   13.33      653.73      770.64      860.49
   0.6283   10.00      598.77      857.80     1046.13
   0.7854    8.00      497.47      582.85      550.77
   0.9425    6.67      367.72      266.33      186.98
   1.0996    5.71      240.65      121.46      104.79
   1.2566    5.00      142.41       76.17       76.74
   1.4137    4.44       81.28       54.20       47.19
   1.5708    4.00       49.13       40.16       41.39
   1.7279    3.64       32.57       27.58       26.46
   1.8850    3.33       22.44       16.52       14.40
   2.0420    3.08       15.53       10.93        9.87
   2.1991    2.86       11.19        8.30        8.32
   2.3562    2.67        8.66        6.18        5.86
   2.5133    2.50        6.93        4.75        4.22
   2.6704    2.35        5.51        4.62        4.35
   2.8274    2.22        4.47        4.91        5.24
   2.9845    2.11        3.61        4.23        4.75
   3.1416    2.00        2.62        2.44        2.27

Spectral Density Using the Bartlett-Priestley Window
Frequency  Period      M = 10      M = 20      M = 30
   0.1571   40.00      604.34      712.73      757.61
   0.3142   20.00      564.28      176.81      107.08
   0.4712   13.33      767.63      927.14      981.10
   0.6283   10.00      900.32     1190.30     1172.23
   0.7854    8.00      607.45      494.85      571.65
   0.9425    6.67      237.16      127.65       87.36
   1.0996    5.71      103.34      113.93      135.34
   1.2566    5.00       75.74       74.88       57.57
   1.4137    4.44       52.64       44.98       38.59
   1.5708    4.00       38.50       44.56       50.59
   1.7279    3.64       27.35       25.28       21.76
   1.8850    3.33       15.68       13.84       13.10
   2.0420    3.08       10.33        9.79        7.41
   2.1991    2.86        7.95        8.31        8.67
   2.3562    2.67        6.04        5.86        7.08
   2.5133    2.50        4.56        3.67        2.90
   2.6704    2.35        4.44        4.38        4.06
   2.8274    2.22        4.99        5.62        5.40
   2.9845    2.11        4.31        5.07        5.08
   3.1416    2.00        2.43        2.23        2.44



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