|
Estimates the nonnormalized spectral density of a stationary time series using a spectral window given the time series data.
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.
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).
Generic: CALL SSWD (X, F, M, PM, SM [, ])
Specific: The specific interface names are S_SSWD and D_SSWD.
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.
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.
where FM(θ) corresponds to the Fejιr kernel of order M.
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.
where M is even. If M is odd, then M + 1 is used instead of M in the above formula.
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, , nf
where nf = 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.
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.
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
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
PHONE: 713.784.3131 FAX:713.781.9260 |