TIMESERIES Procedure

Spectral Density Analysis

Spectral analysis can be performed on the working series by specifying the OUTSPECTRA= option or by specifying the PLOTS=PERIODOGRAM or PLOTS=SPECTRUM option in the PROC TIMESERIES statement. PROC TIMESERIES uses the finite Fourier transform to decompose data series into a sum of sine and cosine terms of different amplitudes and wavelengths. The finite Fourier transform decomposition of the series x Subscript t is

StartLayout 1st Row 1st Column x Subscript t 2nd Column equals StartFraction a 0 Over 2 EndFraction plus sigma-summation Underscript k equals 1 Overscript upper K minus 1 Endscripts f Subscript k Baseline left-parenthesis a Subscript k Baseline cosine omega Subscript k Baseline t plus b Subscript k Baseline sine omega Subscript k Baseline t right-parenthesis 2nd Row 1st Column f Subscript k 2nd Column equals StartLayout Enlarged left-brace 1st Row 1st Column 1 slash 2 2nd Column if upper T is even and k equals upper K minus 1 2nd Row 1st Column 1 2nd Column otherwise EndLayout EndLayout

where

t

is the time subscript, t equals 0 comma 1 comma 2 comma ellipsis comma upper T minus 1

x Subscript t

are the equally spaced time series data

T

is the number of observations in the time series

K

is the number of frequencies in the Fourier decomposition: upper K equals StartFraction upper T plus 2 Over 2 EndFraction if T is even, upper K equals StartFraction upper T plus 1 Over 2 EndFraction if T is odd

k

is the frequency subscript, k equals 0 comma 1 comma 2 comma ellipsis comma upper K minus 1

a 0

is the mean term: a 0 equals 2 x overbar

a Subscript k

are the cosine coefficients

b Subscript k

are the sine coefficients

omega Subscript k

are the Fourier frequencies: omega Subscript k Baseline equals StartFraction 2 pi k Over upper T EndFraction

The Fourier decomposition is performed after the ACCUMULATE=, DIF=, SDIF=, and TRANSFORM= options in the ID and VAR statements have been applied.

Functions of the Fourier coefficients a Subscript k and b Subscript k can be plotted against frequency or against wavelength to form periodograms. The amplitude periodogram upper I Subscript k is defined as follows:

upper I Subscript k Baseline equals StartFraction upper T Over 2 EndFraction left-parenthesis a Subscript k Superscript 2 Baseline plus b Subscript k Superscript 2 Baseline right-parenthesis

Since the Fourier transform is an even, periodic function of frequency that repeats every T ordinates, the periodogram is also. Values of upper I Subscript k for all k therefore can be mapped to the unique values upper I Subscript k Baseline colon k equals 0 comma ellipsis comma upper K minus 1 using the equations

StartLayout 1st Row 1st Column upper I Subscript k 2nd Column equals upper I Subscript negative k Baseline 3rd Column Blank 4th Column for all k 2nd Row 1st Column upper I Subscript k 2nd Column equals upper I Subscript k plus n upper T Baseline 3rd Column Blank 4th Column for n equals plus-or-minus 1 comma plus-or-minus 2 comma plus-or-minus 3 comma ellipsis 3rd Row 1st Column upper I Subscript k 2nd Column equals upper I Subscript upper T minus k Baseline 3rd Column Blank 4th Column for 0 less-than-or-equal-to k less-than-or-equal-to upper K minus 1 EndLayout

The periodogram, upper I Subscript k, is an estimate at the discrete frequencies omega Subscript k of the spectral density function which characterizes the series x Subscript t. By smoothing the periodogram an improved spectral density estimate with reduced variance and bias can be achieved at these points. Smoothing can be accomplished either through use of a spectral window smoothing function or by applying a lag window filter to the series autocovariance function.

When the SPECTRA statement’s DOMAIN=FREQUENCY option is in effect spectral density estimates are computed by smoothing the periodogram ordinates using the equation

StartLayout 1st Row  upper S Subscript k Baseline left-parenthesis upper M right-parenthesis equals sigma-summation Underscript tau equals upper K minus upper T Overscript upper K minus 1 Endscripts w left-parenthesis StartFraction tau Over upper M EndFraction right-parenthesis upper I Subscript k plus tau Baseline EndLayout

where normal w left-parenthesis theta right-parenthesis is the spectral window function whose form is specified by either the KERNEL= option or the WEIGHTS option. M is the kernel scale parameter which acts as a frequency scaling factor in the spectral window smoothing function. Values of upper I Subscript k plus tau that fall outside of 0 less-than-or-equal-to k plus tau less-than-or-equal-to upper K minus 1 are mapped to values inside this range by the equations presented previously.

When the DOMAIN=TIME option is specified, spectral density values are estimated by applying a lag window filter, lamda left-parenthesis h comma upper M right-parenthesis, to the series autocovariance function. The spectral density estimate then can be computed from the filtered autocovariance function using the equation

upper S Subscript k Baseline left-parenthesis upper M right-parenthesis equals sigma-summation Underscript h equals minus left-parenthesis upper T minus 1 right-parenthesis Overscript upper T minus 1 Endscripts lamda left-parenthesis h comma upper M right-parenthesis ModifyingAbove gamma With caret left-parenthesis h right-parenthesis cosine h omega Subscript k

In this case the kernel scale parameter, M, serves as a scale factor for the lag length, h, in the time domain. In the lag window formulation the spectral density estimate is a consistent estimator as upper T comma upper M right-arrow normal infinity under the conditions lamda left-parenthesis h comma upper M right-parenthesis equals 0 for StartAbsoluteValue h EndAbsoluteValue greater-than upper M, and limit Underscript upper T right-arrow normal infinity Endscripts upper M slash upper T equals 0. These conditions lead to the following parameterization of M provided by the SPECTRA statement,

upper M equals c upper K Superscript e

where the values c greater-than 0 and 0 less-than e less-than 1 satisfy the consistency conditions. To specify the kernel scale parameter explicitly, set c equals to the desired scale factor and e equals 0.

For uniformity and computational efficiency, all spectral density estimates are calculated using a spectral window weighting function, w left-parenthesis theta right-parenthesis, applied to the periodogram ordinates. In the case where the DOMAIN=TIME option is specified, the effective spectral window weighting function is computed by the equation

w Subscript TIME Baseline left-parenthesis theta right-parenthesis equals sigma-summation Underscript h equals minus left-parenthesis upper T minus 1 right-parenthesis Overscript upper T minus 1 Endscripts lamda left-parenthesis h comma upper M right-parenthesis cosine h theta

Because the kernel scale parameter, M, serves as a lag scale factor in the time domain and bandwidth scale factor in the frequency domain, the impact of M on spectral density estimates depends on the value of the DOMAIN= option. When DOMAIN=FREQUENCY, increasing values of M decrease variance and increase bias in the spectral density estimates; when DOMAIN=TIME, increasing values of M increase variance and decrease bias.

Using Kernel Specifications

You can specify one of ten different kernel smoothing functions in the SPECTRA statement. Five smoothing functions are available as KERNEL= options, and five complementary smoothing functions, which correspond to lag window filters, are available when the KERNEL= option is used in conjunction with the DOMAIN=TIME option.

For example, a Parzen kernel with a support of 11 periodogram ordinates in the frequency domain can be specified using the kernel option:

spectra / parzen c=5 expon=0;

The TIMESERIES procedure supports the following spectral window kernel functions in the frequency domain where x equals tau slash upper M:

BARTLETT: Bartlett kernel

StartLayout 1st Row 1st Column w left-parenthesis x right-parenthesis 2nd Column equals 3rd Column StartLayout Enlarged left-brace 1st Row 1st Column 1 minus StartAbsoluteValue x EndAbsoluteValue 2nd Column StartAbsoluteValue x EndAbsoluteValue less-than-or-equal-to 1 2nd Row 1st Column 0 2nd Column normal o normal t normal h normal e normal r normal w normal i normal s normal e EndLayout EndLayout

PARZEN: Parzen kernel

StartLayout 1st Row 1st Column w left-parenthesis x right-parenthesis 2nd Column equals 3rd Column StartLayout Enlarged left-brace 1st Row 1st Column 1 minus 6 StartAbsoluteValue x EndAbsoluteValue squared plus 6 StartAbsoluteValue x EndAbsoluteValue cubed 2nd Column 0 less-than-or-equal-to StartAbsoluteValue x EndAbsoluteValue less-than-or-equal-to one-half 2nd Row 1st Column 2 left-parenthesis 1 minus StartAbsoluteValue x EndAbsoluteValue right-parenthesis cubed 2nd Column one-half less-than-or-equal-to StartAbsoluteValue x EndAbsoluteValue less-than-or-equal-to 1 3rd Row 1st Column 0 2nd Column normal o normal t normal h normal e normal r normal w normal i normal s normal e EndLayout EndLayout

QS: quadratic spectral kernel

StartLayout 1st Row 1st Column normal w left-parenthesis x right-parenthesis 2nd Column equals 3rd Column StartFraction 3 Over left-parenthesis 2 pi x right-parenthesis squared EndFraction left-parenthesis StartFraction sine 2 pi x Over 2 pi x EndFraction minus cosine 2 pi x right-parenthesis EndLayout

TUKEY: Tukey-Hanning kernel

StartLayout 1st Row 1st Column w left-parenthesis x right-parenthesis 2nd Column equals 3rd Column StartLayout Enlarged left-brace 1st Row 1st Column left-parenthesis 1 plus cosine left-parenthesis pi x right-parenthesis right-parenthesis slash 2 2nd Column StartAbsoluteValue x EndAbsoluteValue less-than-or-equal-to 1 2nd Row 1st Column 0 2nd Column normal o normal t normal h normal e normal r normal w normal i normal s normal e EndLayout EndLayout

TRUNCAT: truncated kernel

StartLayout 1st Row 1st Column w left-parenthesis x right-parenthesis 2nd Column equals 3rd Column StartLayout Enlarged left-brace 1st Row 1st Column 1 2nd Column StartAbsoluteValue x EndAbsoluteValue less-than-or-equal-to 1 2nd Row 1st Column 0 2nd Column normal o normal t normal h normal e normal r normal w normal i normal s normal e EndLayout EndLayout

When the DOMAIN=TIME option is specified the five kernel functions above are interpreted as lag window filters on the autocovariance function. The lag window kernel functions correspond to the following spectral window smoothing functions where theta equals 2 pi tau slash upper T:

BARTLETT: Bartlett equivalent lag window filter

StartLayout 1st Row 1st Column w left-parenthesis theta right-parenthesis 2nd Column equals 3rd Column StartFraction 1 Over 2 pi upper M EndFraction left-parenthesis StartFraction sine left-parenthesis upper M theta slash 2 right-parenthesis Over sine left-parenthesis theta slash 2 right-parenthesis EndFraction right-parenthesis squared EndLayout

PARZEN: Parzen equivalent lag window filter

StartLayout 1st Row 1st Column w left-parenthesis theta right-parenthesis 2nd Column equals 3rd Column StartFraction 6 Over pi upper M cubed EndFraction left-parenthesis StartFraction sine left-parenthesis upper M theta slash 4 right-parenthesis Over sine left-parenthesis theta slash 2 right-parenthesis EndFraction right-parenthesis Superscript 4 Baseline left-parenthesis 1 minus two-thirds sine squared left-parenthesis theta slash 2 right-parenthesis right-parenthesis EndLayout

QS: quadratic spectral equivalent lag window filter

StartLayout 1st Row 1st Column w left-parenthesis theta right-parenthesis 2nd Column equals 3rd Column StartLayout Enlarged left-brace 1st Row 1st Column StartFraction 3 upper M Over 4 pi EndFraction left-parenthesis 1 minus left-parenthesis upper M theta slash pi right-parenthesis squared right-parenthesis 2nd Column StartAbsoluteValue theta EndAbsoluteValue less-than-or-equal-to pi slash upper M 2nd Row 1st Column 0 2nd Column StartAbsoluteValue theta EndAbsoluteValue greater-than pi slash upper M EndLayout EndLayout

TUKEY: Tukey-Hanning equivalent lag window filter

StartLayout 1st Row 1st Column normal w left-parenthesis theta right-parenthesis 2nd Column equals 3rd Column one-fourth upper D Subscript upper M Baseline left-parenthesis theta minus pi slash upper M right-parenthesis plus one-half upper D Subscript upper M Baseline left-parenthesis theta right-parenthesis plus one-fourth upper D Subscript upper M Baseline left-parenthesis theta plus pi slash upper M right-parenthesis 2nd Row 1st Column upper D Subscript upper M Baseline left-parenthesis theta right-parenthesis 2nd Column equals 3rd Column StartFraction 1 Over 2 pi EndFraction StartFraction sine left-bracket left-parenthesis upper M plus 1 slash 2 right-parenthesis theta right-bracket Over sine left-parenthesis theta slash 2 right-parenthesis EndFraction EndLayout

TRUNC: truncated equivalent lag window filter

StartLayout 1st Row 1st Column normal w left-parenthesis theta right-parenthesis 2nd Column equals 3rd Column upper D Subscript upper M Baseline left-parenthesis theta right-parenthesis EndLayout

Using Specification of Weight Constants

Any number of weighting constants can be specified. The constants are interpreted symmetrically about the middle weight. The middle constant (or the constant to the right of the middle if an even number of weight constants is specified) is the relative weight of the current periodogram ordinate. The constant immediately following the middle one is the relative weight of the next periodogram ordinate, and so on. The actual weights used in the smoothing process are the weights specified in the WEIGHTS option scaled so that they sum to 1.

The moving average calculation reflects at each end of the periodogram to accommodate the periodicity of the periodogram function.

For example, a simple triangular weighting can be specified using the following WEIGHTS option:

spectra / weights 1 2 3 2 1;

Computational Method

If the number of observations, T, factors into prime integers that are less than or equal to 23, and the product of the square-free factors of T is less than 210, then the procedure uses the fast Fourier transform developed by Cooley and Tukey (1965) and implemented by Singleton (1969). If T cannot be factored in this way, then the procedure uses a Chirp-Z algorithm similar to that proposed by Monro and Branch (1977).

Missing Values

Missing values are replaced with an estimate of the mean to perform spectral analyses. This treatment of a series with missing values is consistent with the approach used by Priestley (1981).

Last updated: June 19, 2025