SPECTRUM—Compute and graph spectral estimates

[TomDoan]

Computes and optionally graphs the estimated spectrum of a series.

spectrum.src

Detailed description

Example

Code: Select all

*
* Brockwell & Davis, "Time Series Analysis: Theory and Methods"
* Example 10.4.3 on page 354
*
open data sunspots.dat
calendar 1770
data(format=free,org=columns) 1770:1 1869:1 sunspots
*
diff(center) sunspots / cspots
@spectrum(footer="Figure 10.7 Periodogram of Wolfer Sunspot Numbers",$
   periodogram,nologscale) cspots
@spectrum(footer="Figure 10.8 Spectral Estimate of Wolfer Sunspot Numbers",$
   weights=||3.0,3.0,2.0,1.0||,nologscale,spectrum=smoothed,ordinates=128) cspots
*
* Confidence bands (text only does these for the artificial data). %EDF
* is the Equivalent Degrees of Freedom of the window used by @spectrum.
*
set lower = smoothed*%edf/%invchisqr(.025,%edf)
set upper = smoothed*%edf/%invchisqr(.975,%edf)
graph(footer="Spectral Estimate of Wolfer Sunspot Numbers with 95% Confidence Interval") 3
# smoothed
# lower / 2
# upper / 2
*
* Same thing done on log scale - the confidence bands are equal width
*
graph(footer="Log Spectral Estimate of Wolfer Sunspot Numbers with 95% Confidence Interval",log=10) 3
# smoothed
# lower / 2
# upper / 2

sunspots.dat

Data file for example

(393 Bytes) Downloaded 1163 times

[Future]

Hi,Tom:
I'm a little confused.
In the SPECTRUM.src, @SPECTRUM options, SPECTRUM=(output) series for [0,pi] spectrum (0 frequency at entry 1) .
but in your example,@spectrum(footer="Figure 10.8 Spectral Estimate of Wolfer Sunspot Numbers",$
weights=||3.0,3.0,2.0,1.0||,nologscale,spectrum=smoothed,ordinates=128) cspots.
could you detailed comment like other options?
think you!

[TomDoan]

Hi,Tom:
I'm a little confused.
In the SPECTRUM.src, @SPECTRUM options, SPECTRUM=(output) series for [0,pi] spectrum (0 frequency at entry 1) .
but in your example,@spectrum(footer="Figure 10.8 Spectral Estimate of Wolfer Sunspot Numbers",$
weights=||3.0,3.0,2.0,1.0||,nologscale,spectrum=smoothed,ordinates=128) cspots.
could you detailed comment like other options?
think you!

I'm not sure what your question is. @SPECTRUM is the name of the procedure. It has a SPECTRUM option which saves a series for the estimated spectrum in case you want to do something else with it. In this case it's saved into SMOOTHED which is used to generate upper and lower bounds and do a new graph that includes those.

[kdwest]

I have a question about “spectrum.src”, a procedure I got from the Estima web site.

When using this procedure, I get different answers if I invoke it to (1)compute the periodogram, or (2)“smooth” the periodogram with no smoothing:

(1) @spectrum(logscale=0,periodogram=1) series start end

(2) @spectrum(logscale=0,periodogram=0,width=1,window=flat) series start end

But the two should be the same, according to standard time series and the formula for S(j) given in the description of the WINDOW instruction.

My real question is about how RATS spectral software works, not how this procedure works. I used spectrum.src just to check my understanding.

Can you clarify why the two give different results?

Thanks.

[TomDoan]

With PERIODOGRAM, the procedure doesn't pad the ordinates so they will keep their (theoretical) independence. With any type of window (even width one), they get padded to a useful number (typically a power of 2 times the periods per year so the exact seasonal frequencies are included). You can override the default behavior with an explicit ORDINATES option.

[kdwest]

Thanks. Dumb question. What is the purpose of padding?

[TomDoan]

Historically, speed. The straight FT is an O(T^2) operation, while the Fast Fourier Transform (padding to power of 2) is O(T log2 T). For T=4095 vs T=4096, that's a factor of over 300, which in the 1960's would mean 2 seconds vs 10 minutes. Now it's microseconds vs milliseconds which won't matter much unless you have to do millions of FT's. Now (as I mentioned) we pad by default to the next power of 2 multiplied by the data frequency, so (for instance) 500 odd monthly observations would pad to 12 x 64 = 768. That makes all the harmonics of the seasonal exactly represented as one of the ordinates.