Periodograms

The periodogram of a series is its squared Fourier transform. The cross-periodogram of two series is the product of one FT by the other’s conjugate. You compute both of these using FFT followed by CMULTIPLY.

 

fft 1

fft 2

cmultiply(scale=1./(2*%pi*points)) 1 1 / 3

cmultiply(scale=1./(2*%pi*points)) 1 2 / 4

 

The correct scale factor is \({1 \mathord{\left/{\vphantom {1 {\left( {2{\kern 1pt} \pi {\kern 1pt} N} \right)}}} \right.} {\left( {2{\kern 1pt} \pi {\kern 1pt} N} \right)}}\) where N is the number of actual data points. We represent \(N\) above as POINTS. This might be a number or a variable previously computed.

 

The raw periodogram is used directly in the cumulated periodogram test (see CACCUMULATE), and can be inverse transformed to get the covariogram (see description of lag windows below). However, usually it is smoothed to get a consistent estimate of the spectral density or cross-spectral density.

Smoothing

While the periodogram gives an unbiased estimate of the spectrum, it is inconsistent, as the variance of the estimators does not go to zero as the number of data points grows.

 

RATS provides the method of spectral window smoothing for producing consistent estimates. The alternative procedure is the lag window or weighted covariance method.

 

The instruction WINDOW carries out the smoothing process. This estimates the spectral (cross-spectral) density at an ordinate by taking a weighted average of the periodogram at neighboring ordinates. For a continuous spectral density, there are two requirements for this to produce consistency:

•The window width must go to infinity as the number of data points increases (ensuring that the variance goes to zero). In a practical sense, this means that a window which is too narrow will produce imprecise estimates.

•The window width must increase at a rate slower than the increase in the number of data points (ensuring that the bias goes to zero).

Translation: a very wide window will flatten the peaks and troughs too much.

 

Continuing from above, the following smooths 3 and 4 with a flat window of width 13, producing series 5 and 6.

 

window(width=13) 3 / 5

window(width=13) 4 / 6

 

Tapers

The smoothing method, particularly with the default flat window, can produce some spurious ripples in the estimate due to “data window leakage.” The data set is, theoretically, just a piece of an infinite data set, and we cut it off sharply at both ends. Data window leakage is the result of Fourier transforming this sharp-edged window.

 

A taper reduces this by scaling the ends of the data so they merge smoothly with the zeros on either side. Apply the RATS instruction TAPER to the unpadded part of the series prior to the FFT instruction. Tapering has little effect when you use WINDOW with TYPE=TENT. When you apply a taper, you need to change the scale factor on CMULTIPLY to 1.0/(2*%PI*%SCALETAP), where %SCALETAP is computed by TAPER.

 

Lag Windows

Lag window or weighted covariance estimators were the favored method for computing spectral densities before the Fast Fourier Transform was developed. They are described in Koopmans (1974) on pages 325-6. A simple form of this is used in modern econometrics in the Newey-West covariance estimator—the triangular weighting scheme used in Newey-West is known as the Bartlett window in spectral analysis.

 

The periodogram and the covariogram are inverses of each other under Fourier transformation. The periodogram gives inconsistent estimates of the spectrum because it weights equally (in Fourier transformation) the covariances with a small number of lags, which are estimated with almost a full set of data, and those at high lags, which are computed with almost none. In lag window estimation, we drop a window over the covariogram which gives zero weight to the high lags. The requirement for consistency here is that the number of covariances given non-zero weight goes to infinity at a rate slower than the increase of the number of data points.

 

If you want to employ the lag window method, the simplest way to compute the covariogram is to inverse transform the periodogram. You must pad to at least double the number of data points. The “windowing” process is done by setting up a weight series with the desired form. Make sure to symmetrize it by making the end terms match up with the starting ones. This does a Bartlett window of 13 lags applied to a data set with 130 data points (padded to 320). Series 2 will be the estimate of the spectral density. We sneak the \(1/2\pi\) factor in when we multiply 1 by 2.

 

compute k=13.0

frequency 2 320

rtoc

#  X

#  1

fft 1                               These three lines transform

cmult(scale=1./130) 1 1             series 1 into the covariogram 

ift 1                               of the input series.

cset 2 = %if(t<=k,1-(t-1)/k,0.0)

@cseriessymm 2

cset 2 = (1.0/(2*%pi))*%z(t,1)*%z(t,2)

fft 2

 

Cross-Spectral Statistics

The examples at the start of the section constructed an estimate of the cross-spectral density. Usually, we are interested not in the density itself (which is complex-valued) but statistics derived from it:

 

          

RATS provides the instruction POLAR for computing these statistics. It does a polar decomposition of a series. The phase can be obtained directly from the cross-spectral density, but the coherence requires that we normalize the density by the product of the square roots of the spectral densities of the series. Continuing the example:

 

cmultiply(scale=1./(2*%pi*points)) 2 2 / 7             need spectrum of 2

window(width=13) 7 / 8

cset 7 = %z(t,7)/%csqrt(%z(t,6)*%z(t,8))

polar(periods) 7 / 8 9

 

Series 8 is the coherence, 9 is the phase lead of series 1 over 2 (since we CMULTIPLYed in that order). The PERIODS option converts the phase lead from radians to periods. Note that the phase lead is almost meaningless (and very poorly estimated) when the coherence is small.

 

Procedures

The procedures @SPECTRUM and @CROSSPEC compute and (optionally) graph a spectral density, and coherence and phase, respectively.

 

@crosspec( options )  series1  series2   start  end

@spectrum( options )  series   start  end

 

The series input to the procedures are REAL series. Most of the options are common to both:

 

ORDINATES=Number of ordinates [depends upon length,seasonal]

 

WINDOW=[FLAT]/TENT/QUADRATIC

WIDTH=window width [depends upon ordinates]

 

TAPER=TRAPEZOIDAL/COSINE/[NONE]

WTAPER=Taper width [.25 of length]

 

[GRAPH]/NOGRAPH

HEADER="Header for graph"

FOOTER="Footer for graph"

 

The following segment is taken from the example file SPECTRUM.RPF. The first of these computes the periodogram (unsmoothed) estimate, the second one uses a custom smoothing window of width 7—that requires 4 (relative) weight values for the center and three terms out.

 

@spectrum(footer="Periodogram of Wolfer Sunspot Numbers",$

   periodogram,nologscale) cspots

@spectrum(footer="Spectral Estimate of Wolfer Sunspot Numbers",$

   weights=||3.0,3.0,2.0,1.0||,nologscale,spectrum=smoothed,$

   ordinates=128) cspots

 

 

@ARMASpectrum produces a graph of the spectral density for an input ARMA model where the model is in the form of an EQUATION. This could be an EQUATION defined by BOXJENK or LINREG, or it could be a test equation created using EQUATION:

 

@ARMASpectrum( options )   equation

 

with the options

 

ORDINATES=# of frequency ordinates to use  [512]

HEADER=Header for high-resolution graph

FOOTER=Footer for high-resolution graph

 

 

The @SSMSpectrum procedure returns the estimated multivariate spectrum from a state space model given by the A and SW options. It returns in its argument a series of complex matrices.

 

@SSMSpectrum( options )   spectrum

 

where spectrum will be a SERIES[CMATRIX] created by the procedure with the estimated spectral density. 0 frequency will be in entry 1; in general, frequency \({{2\pi \left( {t - 1} \right)} \mathord{\left/{\vphantom {{2\pi \left( {t - 1} \right)} {T}}} \right.} {T}}\) will be in entry \(t\). In order to avoid problems with non-stationary models, the 0 frequency isn’t calculated. The options are

 

A=transition matrix for the state space model

F=loadings from shocks to the states [identity]

SW=covariance matrix of the shocks in the state equation

These are the same options as would be used in a DLM instruction.

 

ORDINATES=# of frequency ordinates used for calculations [512]

This will return values only out to frequency \(\pi\), thus to entries \({{T} \mathord{\left/{\vphantom {{T} 2}} \right.} 2} + 1\) where \(T\) is the choice for ORDINATES.

 

COMPONENTS=VECT[INT] of components 

Use this to indicate which state components (numbered in 1,...,# of states) that you actually want. For instance, COMPONENTS=||1|| will return the estimate just for the first component.

 

Copyright © 2025 Thomas A. Doan