38.0 SPECTRAL Command
The SPECTRAL command computes spectral and cross-spectral
densities of a multivariate time series. The periodogram and cross
periodogram can be calculated. The user can optionally output the cosine
and sin coefficients for each series, the periodogram, the real part of
the cross-periodogram, the imaginary part of the cross-periodogram, the
spectral density estimate, the cospectral density estimate, the
quadrature-spectrum estimate, the amplitude, the coherency squared and
the phase spectrum. The B34S SPECTRAL command calculates many of the
same values as the SAS/ETS(r) PROC SPECTRA. The B34S SPECTRAL command
uses the FORTRAN routines EZFFTI and EZFFTF which were built by Paul
Swarztrauber and are available in the FFTPACK library. The SPECTRAL
command has added capability from the SPECRTAL sentence in the BJIDEN
and BJEST command.
Form of Command.
B34SEXEC SPECTRAL options parameters $
VAR Xvar1 Xvar2 $
WEIGHTS( ) $
B34SEEND$
The VAR sentence is required. Unless the LIST, PLOT or OUTPUT parameters
on the SPECTRAL sentence are supplied, nothing will be done.
SPECTRAL command (sentence) parameters.
IBEGIN=n1 Sets beginning observation. Defaults to 1.
IEND=n2 Sets ending observation. Defaults to NOOB.
PLOT(k1, k2) Sets spectral values to plot using line printer
graphics. Key words for k1 ... are:
SIN, COS, P, S, RP, IP, CS, QS, A, K, PH, ALL.
These key words are defined below.
LIST(k1, k2) Sets spectral values to list using line printer.
Key words for k1 ... are:
SIN, COS, P, S, RP, IP, CS, QS, A, K, PH, ALL.
These key words are defined below.
PLOTBY(k3) Sets whether plots are by FREQ or PERIOD. Valid
values for k3 are PERIOD or FREQ. The default is
FREQ.
SCAFNAME(key) Sets SCA FSAVE DSN name. Up to 8 letters can be
supplied. Default = SPEC.
SCAUNIT=n3 Sets SCA FSAVE file unit. Default = 44.
OUTPUT(k1, k2) Sets spectral values to output to SCA FSAVE file.
Key words for k1 ... are:
SIN, COS, P, S, RP, IP, CS, QS, A, K, PH, ALL.
These key words are defined below. If the OUTPUT
parameter is used, the file for the SCA FSAVE must
have been opened before the SPECTRAL command is
encountered. The SPECTRAL command will add a SCA
FSAVE file with DSN name supplied by SCAFNAME
to the files already on the unit.
NOBSPP=n4 Sets number of observations per period. Default =1.
If this parameter is set, FREQ will be expressed
cycles per unit of time. If NOBSPP=0, FREQ values
will agree with SAS.
VAR sentence.
The VAR sentence specifies the variables to pass to the SPECTRAL
command. If more than two series are passed, CS, QS, A, K and PH are
with respect to the last series. A maximum of 7 series can be listed.
The form of the VAR sentence is:
VAR Xvar1 Xvar2 $
WEIGHTS sentence.
The WEIGHTS sentence allows the user to specify an odd number of
weights to be used for smoothing the spectrum. The weights are
normalized to sum to 1/(4*pi). Up to 99 weights can be supplied. The
number of weights supplied must be odd. The form of the WEIGHTS
sentence is:
WEIGHTS(r1, r2) $
where r1, ... are the weight values
Note: Unless the LIST or PLOT parameter is passed, the SPECRAL command
will not generate any output.
References.
Brocklebank, John and David Dickey. "SAS System for Forecasting
Time Series," Cary, NC: SAS Institute Inc, 1986, 241p.
Hamilton, James "Time Series Analysis," Princeton University
Press: Princeton, New Jersey, 1994.
Jenkins, Gwilym and Donald Watts. "Spectral Analysis and it's
Applications," Holden Day 1968, 524p.
SAS/ETS(r) User's Guide, Version 6, First Edition, Cary, NC: SAS
Institute Inc, 1988, 560p.
Discussion of quantities calculated.
Define N as number of observations. If N = even, K = (N/2) - 1.
If N = odd, K = (N-1)/2. If the WEIGHTS sentence is found,
the supplied p weights are normalized to sum to 1/4*pi and are saved
in w(1) .... w(p) where p is odd.
FREQ - Frequency from 0 to pi if NOBSPP is set = 0. Cycles per
observation is NOBSPP*FREQ/2pi where NOBSPP is the number.
of observations per unit time. Using the default setting of
NOBSPP = 1, FREQ is in the range 0 to .5.
PERIOD - Period or wavelength ( 1 / FREQ ).
COS_n - Cosine transform of Xn. Defined over range 1 to K. COS_n(k)=
sum for i=1,N of (2/N)*Xn(i)*cos(k*(i-1)*2*pi/N).
SIN_n - Sine transform of Xn. Defined over range 1 to K. SIN_n(k)=
sum for i=1,N of (2/N)*Xn(i)*sin(k*(i-1)*2*pi/N).
P_n - Periodogram of Xn. Defined over range 1 to K. P_n(k) =
(N/2)*((COS_n**2) + (SIN_n**2)).
S_n - Spectral density estimate of Xn. Defined over range 1 to K as
the weighted periodogram. S_n(k) = sum for i=-(p-1)/2 to
(p-1)/2 of w(i)*P_n(k+i). If the WEIGHTS sentence is not
present, S_n = P_n.
Note: The following values are only defined if more than one series
is supplied in the VAR sentence.
RP_n_m - Real part of cross periodogram of Xn and Xm. Defined over
range 1 to K as RP_n_m = (N/2)*((COS_n*COS_m)+(SIN_n*SIN_m)).
IP_n_m - Imaginary part of cross periodogram of Xn and Xm. Defined over
range 1 to K as IP_n_m = (N/2)*((COS_n*SIN_m)-(SIN_n*COS_m)).
CS_n_m - Cospectral density estimate (real part of cross spectrum) or
the weighted real part of the cross periodogram. Defined over
range 1 to K as CS_n_m(k) = sum for i=-(p-1)/2 to (p-1)/2 of
w(i)*RP_n_m(k). If the WEIGHTS sentence is not present,
CS_n_m = RP_n_m.
QS_n_m - Quadrature-spectrum (imaginary part of cross spectrum) or
the weighted imaginary part of the cross periodogram. Defined
over range 1 to K as QS_n_m(k) = sum for i=-(p-1)/2 to (p-1)/2
of w(i)*IP_n_m(k). If the WEIGHTS sentence is not present,
QS_n_m = IP_n_m.
A_n_m - Amplitude (modulus of cross-spectrum). Defined over range 1
to K as A_n_m(k) = ((CS_n_m**2) + (QS_n_m**2))**.5.
K_n_m - Coherency squared. Defined over range 1 to K as
K_n_m(k)=(A_n_m**2)/(S_n*S_m). If the WEIGHTS sentence is
not present, K_n_m reduces to 1.0 for all frequencies.
PH_n_m - Phase spectrum in radiams. Defined over range 1 to K as
PH_n_m(k) = arctan(QS_n_m , CS_n_m)
Note: To simplify notation and assuming k series were loaded, B34S
defines RP_i_k as RP_i, IP_i_k as IP_i, CS_i_k as CS_i,
QS_i_k as QS_i, A_i_k as A_i, K_i_k as k_i and PH_i_k as PH_i.
Note: For smoothing at the end points it is assumed that the moving
average reflects and the end point periodogram values are reused
as many times as needed.
If WEIGHT sentence is not used, K_n_m = 1.0 at all frequencies.
Example 1.
Assuming the Gas Data has been loaded, the following example calculates
the sin and cosine transforms, the periodogram and the spectrum. All
values are saved in the SCA FSAVE file spec.fsv in dsn gasspec.
b34sexec options open('spec.fsv') unit(44) disp=unknown$ b34srun$
b34sexec options clean(44)$ b34srun$
b34sexec spectral scafname='gasspec' scaunit=44
output(sin cos p s)$
weights(1 1 1)$
var gasin $
b34srun$
Example 2.
Assuming the Gas Data has been loaded, the following example calculates
the periodogram and the spectrum for each series and the amplitude,
coherency squared and the phase spectrum. Triangular weights are used.
All values are saved in SCA FSAVE file SPEC.FSV in dsn GASSPEC.
b34sexec options open('spec.fsv') unit(44) disp=unknown$ b34srun$
b34sexec options clean(44)$ b34srun$
b34sexec spectral scafname='gasspec' scaunit=44
output(p s a k ph)$
weights(1 2 1)$
var gasin gasout$
b34srun$
Example 3.
The periodogram can be calculated with OLS. The below listed code
uses the B34S MACRO olsspec to calculate the periodogram using
the B34S command RR. These values are tested against those obtained
with the B34S command SPECTRAL. This example replicates Brocklebank
and Dickey (1986) page 201. Odd and even numbers of observations
are tested.
b34sexec options include('c:\b34slm\b34sdata.mac') macro(water)$
b34srun$
b34sexec spectral list(sin,cos,p,s) nobspp=0$ var Y$ b34srun$
%b34smacro olsspec$
/$ Needs to be called as %b34smcall(y=yname noob=_ ntest=_ )
/$ y = yname
/$ noob = number of observations i series
/$ ntest must be
%b34sdo i=1,&ntest$
b34sexec data set $
%b34sif(&i.eq.1)%then$
build sin cos ii freq$
%b34sendif$
gen freq=(timespi(2.0)/%b34seval(&noob))*%b34seval(&I)$
gen ii =timespi(2.0)/freq$
gen sin=sin(timespi(2.0)*kount()/ii)$
gen cos=cos(timespi(2.0)*kount()/ii)$
b34srun$
b34sexec rr$ model %b34seval(&y)=sin cos$ b34srun$
%b34senddo $
%b34smend$
%b34smcall olsspec(y=y noob=56 ntest=20)$
B34SEXEC OPTIONS INCLUDE('C:\b34slm\b34sdata.mac') macro(WATER)$
B34SRUN$
b34sexec data set$
gen if(kount().eq.56)y=delobs()$
b34seend$
b34sexec spectral list(sin,cos,p,s) nobspp=0$ var Y$ b34srun$
%b34smcall olsspec(y=y noob=55 ntest=20)$