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
