Chapter 8: Time Series Analysis and Forecasting

SWEP

Estimations of the nonnormalized spectral density of a stationary time series based on specified periodogram weights given the periodogram.

Required Arguments

N — Number of observations in the appropriately 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 ωk = 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)

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.

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

NWT — Number of weights.   (Input)
NWT must be greater than or equal to one.
Default: NWT = size (WT,1).

FORTRAN 90 Interface

Generic:                              CALL SWEP (N, PX, F, WT, SX [,…])

Specific:                             The specific interface names are S_SWEP and D_SWEP.

FORTRAN 77 Interface

Single:                                CALL SWEP (N, PX, NF, F, NWT, WT, SX)

Double:                              The double precision name is DSWEP.

Description

Routine SWEP estimates the nonnormalized spectral density function of a stationary time series using a fixed sequence of weights 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

(Here, a means the greatest integer less than or equal to a.) The symmetry of the periodogram is used to recover the ordinates at negative frequencies.

Consider the sequence of m = NWT weights {wj} for j = m/2, …, (m m/2 1) where Σjwj = 1. These weights are fixed in the sense that they do not depend on the frequency ω at which to estimate the nonnormalized spectral density hX(ω). The estimate of the nonnormalized spectral density is computed according to

where

and k(ω) is the integer such that ωk,0 is closest to ω. The weights specified by argument WT may be relative since they are normalized to sum to one in the actual computation of

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. Note that periodogram ordinate

is replaced by

and the sum reflects at each end.

The nonnormalized spectral density estimate is computed over the set of frequencies

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

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.

Approximate confidence intervals for h(ω) can be computed using formulas given in the introduction.

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) 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 SWEP to these data produces the following results:

 

      USE IMSL_LIBRARIES

 

      IMPLICIT   NONE

      INTEGER    LDPM, LDRDAT, N, NDRDAT, NF, NOBS, NPAD, NWT

      PARAMETER  (LDRDAT=176, NDRDAT=2, NF=20, NOBS=100, NWT=7, &

                NPAD=NOBS-1, LDPM=(NOBS+NPAD)/2+1, N=NOBS+NPAD)

!

      INTEGER    I, IFSCAL, IPVER, NROW, NVAR

      REAL       F(NF), PI, PM(LDPM,5), RDATA(LDRDAT,NDRDAT), &

                 FLOAT, SM(NF,2), SX(NF), WT(NWT), X(NOBS)

      CHARACTER  CLABEL(3)*30, FMT*20, RLABEL(1)*4, TITLE*28

      INTRINSIC  FLOAT

!

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

!

      DATA WT/1., 2., 3., 4., 3., 2., 1./

      DATA IPVER/1/, IFSCAL/0/

      DATA FMT/'(F9.4)'/

      DATA CLABEL/'     ', '%/Frequency', 'Spectral%/Estimates'/

      DATA RLABEL/'NONE'/

      DATA TITLE/'Results of Spectral Analysis'/

!                                 Initialization

      PI = 2.0*ASIN(1.0)

      DO 10  I=1, NF

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

   10 CONTINUE

!                                 Wolfer Sunspot Data for years

!                                 1770 through 1869

      CALL GDATA (2, RDATA, NROW, NVAR)

!                                 Compute modified periodogram

      CALL PFFT (X, PM, IPVER=IPVER)

!

!                                 Compute spectral density

      CALL SWEP (N, PM(:,3), F, WT, SX)

!

!                                 Print results

!

!                                 Copy the frequencies to the output

!                                 matrix

      CALL SCOPY (NF, F, 1, SM(1:,1), 1)

!                                 Copy the spectral estimates to the

!                                 output matrix

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

!                                 Call printing routine

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

!

      END

Output

 

Results of Spectral Analysis
            Spectral
Frequency  Estimates
   0.1571   710.8386
   0.3142   116.3940
   0.4712   937.1508
   0.6283  1209.8268
   0.7854   538.9236
   0.9425    84.9561
   1.0996   128.0791
   1.2566    55.0304
   1.4137    40.2022
   1.5708    46.4240
   1.7279    21.0053
   1.8850    12.1449
   2.0420     8.8654
   2.1991     7.2589
   2.3562     6.8078
   2.5133     3.3873
   2.6704     3.9504
   2.8274     5.7418
   2.9845     4.4652
   3.1416     4.1216



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