SSWD
Estimates the nonnormalized spectral density of a stationary time series using a spectral window given the time series data.
Required Arguments
X — Vector of length NOBS containing the stationary time series. (Input)
F — Vector of length NF containing the frequencies at which to evaluate the spectral density estimate. (Input)
The units of F correspond to the scale specified by IFSCAL. The elements of F must be in the range (‑π/TINT,π/TINT), inclusive for IFSCAL = 0 and (‑1/(2 * TINT), 1/(2 * TINT)) inclusive for IFSCAL = 1.
M — Vector of length NM containing the values of the spectral window parameters M. (Input)
For the Parzen spectral window (ISWVER = 5), all values of the spectral window parameters M must be even.
PM — (⌊N/2⌋ + 1) by 5 matrix that contains a summarization of the periodogram analysis. (Output)
For k = 0, 1, …, ⌊N/2⌋, the (k + 1)-st element of the j‑th column of PM 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 |
Periodogram ordinate, I(ωk). |
|
4 |
Cosine transformation coefficient, A(ωk). |
|
5 |
Sine transformation coefficient, B(ωk). |
Note N = NOBS + NPAD.
SM — NF by NM + 2 matrix containing a summarization of the spectral analysis. (Output)
The k‑th element of the j‑th column of SM is defined as
|
Col. |
Description |
|
1 |
Frequency, F(k). |
|
2 |
Period, pk where pk = 2π/F(k) for IFSCAL = 0 and pk = 1/F(k) for IFSCAL = 1. If F(k) = 0, pk is set to missing. |
|
3 |
Spectral density estimate at F(k) using the spectral window parameter M(1). |
|
4 |
Spectral density estimate at F(k) using the spectral window parameter M(2). |
|
⋮ |
|
|
NM +2 |
Spectral density estimate at F(k) using the spectral window parameter M(NM) |
where k = 1, …, NF.
Optional Arguments
NOBS — Number of observations in the stationary time series X. (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 transform and sine transform of the centered and padded time series, and the spectral density estimate based on a specified version of a spectral window for a given set of spectral window parameters. |
XCNTR — Constant used to center the time series X. (Input)
Default: XCNTR = the arithmetic mean.
NPAD — Number of zeroes used to pad the centered time series. (Input)
NPAD must be greater than or equal to zero.
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. |
NF — Number of frequencies at which to evaluate the spectral density estimate. (Input)
Default: NF = size (F,1).
TINT — Time interval at which the series is sampled. (Input)
For a discrete parameter process, usually TINT = 1. For a continuous parameter process, TINT > 0. TINT is used to adjust the spectral density estimate.
Default: TINT = 1.0.
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.
NM — Number of spectral window parameters M used to compute the spectral density estimate for a given spectral window version. (Input)
NM must be greater than or equal to one.
Default: NM = size (M,1).
LDPM — Leading dimension of PM exactly as specified in the dimension statement of the calling program. (Input)
LDPM must be greater than or equal to ⌊N/2⌋ + 1.
Default: LDPM = size (PM,1).
LDSM — Leading dimension of SM exactly as specified in the dimension statement of the calling program. (Input)
LDSM must be greater than or equal to NF.
Default: LDSM = size (SM,1).
FORTRAN 90 Interface
Generic: CALL SSWD (X, F, M, PM, SM [, …])
Specific: The specific interface names are S_SSWD and D_SSWD.
FORTRAN 77 Interface
Single: CALL SSWD (NOBS, X, IPRINT, XCNTR, NPAD, IFSCAL, NF, F, TINT, ISWVER, NM, M, PM, LDPM, SM, LDSM)
Double: The double precision name is DSSWD.
Description
Routine SSWD estimates the nonnormalized spectral density function of a stationary time series using a spectral window given a sample of n = NOBS observations {Xt} for t = 1, 2, …, n.
Let
for t = 1, …, N represent the centered and padded data where N = NOBS + NPAD,
and
is determined by
The modified periodogram of
for t = 1, …, N is estimated by
where
and
represent the
cosine and sine transforms, respectively, and K is the scale factor equal to 1/(2πn). Since the periodogram is an even function of the frequency, it is sufficient to estimate the periodogram over the discrete set of nonnegative frequencies
The routine PFFT is used to compute the modified periodogram of
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 Fejér 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
The argument NM specifies the number of window parameters M and corresponds to the number of spectral density estimates to be computed for a given spectral window. The nonnormalized spectral density is estimated over the set of frequencies
ω = fi i = 1, …, nƒ
where nƒ = NF. These frequencies are in the scale specified by the argument IFSCAL but are transformed to the scale of radians per unit time for computational purposes.
The above formula for
assumes the data {Xt} correspond to a realization of a discrete‑parameter stationary process observed consecutively in time. In this case, the observations are equally spaced in time with interval Δt = TINT equivalent to one. However, if the data correspond to a realization of a continuous‑parameter stationary process recorded at equal time intervals, then the estimate of the nonnormalized spectral density must be adjusted for the effect of aliasing. In general, the estimate of hX(ω) is given by
Note that the frequency ω of the desired spectral estimate is assumed to be input in a form already adjusted for the time interval Δt. Approximate confidence intervals for h(ω) can be computed using formulas given in the introduction.
Comments
1. Workspace may be explicitly provided, if desired, by use of S2WD/DS2WD. The reference is:
CALL S2WD (NOBS, X, IPRINT, XCNTR, NPAD, IFSCAL, NF, F, TINT, ISWVER, NM, M, PM, LDPM, SM, LDSM, CX, COEF, WFFTC, CPY)
The additional arguments are as follows:
CX — Complex vector of length N containing the centered and padded time series X. (Output)
COEF — Complex vector of length N containing the Fourier coefficients of the finite Fourier transform of CX. (Output)
Note that COEF(k) is the appropriately scaled Fourier coefficient at frequency ωk, k = 0, 1, …, N ‑ 1.
WFFTC — Vector of length 4N + 15.
CPY — Vector of length 2N.
2. The normalized spectral density estimate is obtained by dividing the nonnormalized spectral density estimate in matrix SM by an estimate of the variance of X.
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 SSWD to these data produces the following results:
USE GDATA_INT
USE SSWD_INT
USE WRRRL_INT
IMPLICIT NONE
INTEGER LDPM, LDSM, NF, NM, NOBS
REAL PI
PARAMETER (NF=20, NM=3, NOBS=100, PI=3.141592654, &
LDPM=NOBS, LDSM=NF)
!
INTEGER I, ISWVER, M(NM), NCOL, NROW
REAL F(NF), PM(LDPM,5), RDATA(176,2), FLOAT, SM(LDSM,5), &
X(NOBS)
CHARACTER CLABEL(6)*9, FMT*20, RLABEL(1)*6, TITLE*60
INTRINSIC REAL
!
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 (Default)
! USE Default Frequency in radians per unit time
! Determine frequencies at which
! to evaluate spectral density
DO 10 I=1, NF
F(I) = PI*FLOAT(I)/FLOAT(NF)
10 CONTINUE
! USE Default Time interval for discrete data
! Spectral window parameters
M(1) = 10
M(2) = 20
M(3) = 30
! Compute spectral density using
! the Parzen window
ISWVER = 5
CALL SSWD (X, F, M, PM, SM, ISWVER=ISWVER)
! 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
CALL SSWD (X, F, M, PM, SM, ISWVER=ISWVER)
! 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
