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

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