SWED
Estimations of the nonnormalized spectral density of a stationary time series based on specified periodogram weights 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.
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.
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 the missing value or NaN (not a number). |
|
3 |
Periodogram ordinate, I(ωk). |
|
4 |
Cosine transformation coefficient, A(ωk). |
|
5 |
Sine transformation coefficient, B(ωk). |
SM — NF by 3 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 specified weights WT. |
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 and sine transforms of the centered and padded time series, and the spectral density estimate based on a specified weight sequence. |
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. The length of the centered and padded time series is N = NOBS + NPAD.
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)
NF must be greater than zero.
Default: NF = size (F,1).
TINT — Time interval at which the series is sampled. (Input)
For a discrete parameter process, usually TINT = 1.0. For a continuous parameter process, TINT > 0.0. TINT is used to adjust the spectral density estimate.
Default: TINT = 1.0.
NWT — Number of weights. (Input)
NWT must be greater than or equal to one.
Default: NWT = size (WT,1).
LDPM — Leading dimension of PM exactly as specified in the dimension statement in 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 in the calling program. (Input)
LDSM must be greater than or equal to NF.
Default: LDSM = size (SM,1).
FORTRAN 90 Interface
Generic: CALL SWED (X, F, WT, PM, SM [, …])
Specific: The specific interface names are S_SWED and D_SWED.
FORTRAN 77 Interface
Single: CALL SWED (NOBS, X, IPRINT, XCNTR, NPAD, IFSCAL, NF, F, TINT, NWT, WT, PM, LDPM, SM, LDSM)
Double: The double precision name is DSWED.
Description
Routine SWED estimates the nonnormalized spectral density function of a stationary time series using a fixed sequence of weights, given a sample of n = NOBS observations {Xt}, for
t = 1, 2, …, n.
Let
represent the centered and padded data where N = NOBS + NPAD,
and
is determined by
The modified periodogram of
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 at the discrete set of nonnegative frequencies
(Here, ⌊a⌋ means the greatest integer less than or equal to a.) The routine PFFT is used to compute the modified periodogram of
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 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 S2ED/DS2ED. The reference is:
CALL S2ED (NOBS, X, IPRINT, XCNTR, NPAD, IFSCAL, NF, F, TINT, NWT, WT, 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 + 1) is the appropriately scaled Fourier coefficient at frequency ωk, k = 0, 1, …, N ‑ 1.
WFFTC — Work vector of length 4N + 15.
CPY — Work vector of length 2N.
2. The centered and padded time series is defined by
CX(j) = X(j) ‑ XCNTR for j = 1, …, NOBS
CX(j) = 0 for j = NOBS + 1, …, N
where N = NOBS + NPAD.
3. 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 SWED to these data produces the following results:
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER LDPM, LDRDAT, LDSM, NDRDAT, NF, NOBS, NPAD, NWT
PARAMETER (LDRDAT=176, NDRDAT=2, NF=20, NOBS=100, NWT=7, &
LDSM=NF, NPAD=NOBS-1, LDPM=(NOBS+NPAD)/2+1)
!
INTEGER I, NROW, NVAR
REAL F(NF), PI, PM(LDPM,5), RDATA(LDRDAT,NDRDAT), &
REAL, SM(LDSM,3), WT(NWT), X(NOBS), IFSCAL, IPRINT, TINT
CHARACTER CLABEL(4)*20, FMT*20, RLABEL(1)*4, TITLE*28
INTRINSIC FLOAT
!
EQUIVALENCE (X(1), RDATA(22,2))
!
DATA WT/1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0/
DATA IPRINT/0/, IFSCAL/0/, TINT/1.0/
DATA FMT/'(F9.4, F6.2, F9.4)'/
DATA RLABEL/'NONE'/
DATA CLABEL/' ', '%/Frequency', '%/Period', 'Spectral%/Estimates' &
/
DATA TITLE/'Results of Spectral Analysis'/
! Initializations
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)
! Center on arithmetic mean
! Spectral density
CALL SWED (X, F, WT, PM, SM)
! Print Results
CALL WRRRL (TITLE, SM, RLABEL, CLABEL, FMT=FMT)
!
END
Output
Results of Spectral Analysis
Spectral
Frequency Period Estimates
0.1571 40.00 710.8386
0.3142 20.00 116.3940
0.4712 13.33 937.1508
0.6283 10.00 1209.8268
0.7854 8.00 538.9236
0.9425 6.67 84.9561
1.0996 5.71 128.0791
1.2566 5.00 55.0304
1.4137 4.44 40.2022
1.5708 4.00 46.4240
1.7279 3.64 21.0053
1.8850 3.33 12.1449
2.0420 3.08 8.8654
2.1991 2.86 7.2589
2.3562 2.67 6.8078
2.5133 2.50 3.3873
2.6704 2.35 3.9504
2.8274 2.22 5.7418
2.9845 2.11 4.4652
3.1416 2.00 4.1216
