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)$